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ON  THE  CALCULATION  OF  CRITICAL  POINTS  BY  THE  METHOD  OF 

HEIDEMANN  AND  KHALIL 

Brian  E.  Eaton 


The  formulation  of  critical  point  criteria  by  Heidemann  and  Khalil  is  analyzed 
and  contrasted  with  the  original  formulation  of  Gibbs.  An  extension  to  the  solution 
technique  originally  used  by  Heidemann  and  Khalil,  and  later  improved  by  Michelsen 
and  Heidemann,  is  presented  along  with  its  detailed  implementation  for  a  general 
two-constant  cubic  equation  of  state.   Finally,  FORTRAN  software  developed  for  these 
computations  is  carefully  documented. 

Key  words:   algorithm;  critical  locus;  Helmholtz  free  energy;  mixtures;  Peng-Robinson 
equation  of  state;  thermodynamics 

1  .   Introduction 

The  topic  of  this  technical  note  is  the  computation  of  critical  points  in  multicomponent 
mixtures  using  equations  of  state.  The  thermodynamic  conditions  for  a  critical  point  were 
orglnally  formulated  by  Gibbs  \3}i  but  because  of  their  complexity,  they  were  not  solved 
using  equations  of  state  for  multicomponent  mixtures  for  another  100  years.  Peng  and 
Robinson  \10}  computed  critical  points  in  multicomponent  mixtures  using  their  cubic  equation 
of  state  and  the  conditions  for  criticality  exactly  as  given  by  Gibbs  \3},  i.e.,  Gibbs' 
equations  (206-208).   Shortly  thereafter,  Heidemann  and  Khalil  \4}  proposed  an  alternative 
computational  procedure  based  on  a  mathematical  reformulation  of  the  original  conditions  of 
Gibbs.  The  computational  procedure  of  Heidemann  and  Khalil  is  substantially  more  efficient 
than  that  employed  by  Peng  and  Robinson. 

The  purpose  of  this  note  is  twofold.   First,  the  criteria  for  a  critical  point  will  be 
examined  in  some  detail.   Emphasis  is  placed  on  elucidating  the  similarities  and  differences 
between  the  two  mathematical  formulations  mentioned  above.  Second,  a  description  of  Heidemann 
and  Khalil' s  solution  technique  using  a  general  two-constant  cubic  equation  of  state  is 
presented  along  with  an  extension  which  may  lead  to  a  more  efficient  iteration  on  the  volume. 
In  the  final  section,  the  FORTRAN  software  developed  for  the  computations  is  carefully 
documented. 

2.   Conditions  for  the  Local  Stability  of  a  Homogeneous  Phase 

To  formulate  the  mathematical  conditions  for  a  critical  point,  we  begin  by  considering 
the  stability  of  a  homogeneous  phase  with  respect  to  continuous  change.   Implicit  in  this 
statement  is  the  differentiation  between  an  arbitrarily  small  (continuous)  change  in  a  phase 
and  a  change  which  involves  the  appearance  of  a  new  phase  (discontinuous).   A  phase  may  be 
stable  with  regard  to  the  former  type  of  change  and  unstable  with  regard  to  the  latter  (i.e., 
a  metastable  state).   Stability  with  respect  to  continuous  changes  is  referred  to  as  local 
stability  and  that  with  respect  to  discontinuous  changes  is  termed  global  stability. 

Thermodynamic  stability  theory  rests  upon  a  hypothesis  which,  in  its  original  form  as 
stated  by  Gibbs,  does  not  lead  unambiguously  to  intuitive  concepts  and  mathematical  formalism; 
an  excerpt  from  Tisza  \1iJ}  (p.  41)  clarifies  this  point: 

Gibbs  stated  the  extremum  principle  in  two  versions:   in  an  isolated  system  the 
entropy  tends  to  a  maximum  at  constant  energy,  or  alternatively,  the  energy  tends  to 
a  minimum  at  constant  entropy.  Although  these  statements  undoubtedly  express 
important  truths,  they  lack  precision  to  the  point  of  being  paradoxical.  If  an 
isolated  system  is  not  in  equilibrium,  we  can  associate  no  entropy  with  it,  and  if 


it  is  in  equilibrium,  its  entropy  can  no  longer  increase.  Many  authors  have 
grappled  with  this  dilemma  until  a  satisfactory  solution  was  found  in  terms  of  the 
composite  system.  Consider  a  system  of  two  or  more  spatially  disjoint  parts 
separated  by  adiabatic  partitions.   After  reaching  equilibrium  each  part  has  a 
definite  entropy,  the  sum  of  which  is  associated  with  the  composite  system.   The 
relaxation  of  an  internal  constraint  will,  in  general,  trigger  a  process,  namely  the 
redistribution,  say,  of  energy  at  constant  volume, 

(DU-)^,  +  (DU")^„    =   0 

leading  to  a  new  equilibrium  with  higher  entropy.  Here  D  denotes  a  finite  change 
that  actually  takes  place  in  the  system  and  for  which  all  the  conservation  laws  and 
boundary  conditions  are  satisfied.   Thus  the  increase  of  entropy  is  perfectly  well 
defined  since  it  is  associated  with  a  transition  from  a  more  restricted  to  a  less 
restricted  equilibrium. 

Development  of  the  stability  conditions  from  the  extremum  principles  and  the  composite 
system  concept  summarized  below  can  also  be  found  in  Tisza  [14]  and  in  Modell  and  Reid  [8]. 
It  is  convenient  to  consider  the  results  in  terms  of  the  energy  rather  than  the  entropy  since 
the  Legendre  transforms  of  the  energy  result  in  the  free  energy  expressions  which  are  readily 
calculated  using  equations  of  state. 

The  extensive  state  of  a  system  can  be  represented  in  terms  of  c  +  2  variables,  where  c 
is  the  number  of  independent  components, 

U(S,  V,  n^ n^_^  ,  N);  (2.1) 

the  internal  energy  (U)  is  a  function  of  the  entropy  (S),  volume  (V),  number  of  moles  of 
component  i  (n^)  and  the  total  number  of  moles  (N).   It  is  generally  assumed  that  the 
stability  of  a  system  is  an  intrinsic  property  which  does  not  depend  on  its  size.   This 
assumption  is  not  universally  valid:   for  example,  it  is  not  valid  for  small  droplets,  where 
surface  tension  effects  are  important,  nor  is  it  valid  for  large  stars,  where  gravity  must  be 
taken  into  account.   However,  size  dependent  effects  will  not  be  considered  here;  hence,  the 
size  of  the  system  will  be  fixed  by  holding,  say,  N  constant.   Variations  in  the  energy  are 
then  expressed  by 

1=1    1 

where  the  variables  held  constant  in  the  partial  differentiations  follow  from  the  functional 
dependence  expressed  in  (2.1).  Introducing  some  additional  nomenclature,  we  rewrite  (2.2)  as 
follows, 

dU   =    I       f-dX.   =    I        U.  dX.   ,  (2.3) 

1=1     1         1=1 

where  X^  is  the  i^*^  component  of  the  vector  X  =  [S,  V,  nj^,  ...,  ng--]]'^, 

U.      .     —      ,  and     r  .  c  .  1  . 

The  Uj_  represent  the  intensities  of  the  system,  that  Is  U^  is  the  i'''^  component  of  the  vector 
P  =  [T,  -p,  PI,  ...,  Uq-i  ]T  where  T  is  temperature,  p  is  pressure,  and  \ii   is  the  chemical 
potential  of  mixture  component  i.   The  superscript  T  represents  the  transpose  operation,  i.e., 
X  and  P  are  column  vectors. 


The  extremum  principle  is  expressed  as  follows.   Consider  a  composite  system,  the 
subsystems  of  which  are  identified  by  the  superscripts  a   and  6.   The  initial  state,  in  which 
subsystems  a   and  B  are  assumed  to  be  in  a  state  of  unconstrained  equilibrium  (all  the 
intensities  in  a  equal  the  intensities  of  g) ,  is  characterized  by 

^jsystem  ^  ^^a  ^  ^B   _  (2_^) 

The  constraint  that  the  composite  system  be  isolated  takes  the  form 

x"  +  X^  =  const.,     i  =  1,  ....  r  (2.5) 

C(      S 

and  the  relative  sizes  of  a  and  6  are  specified  by  fixing  the  valves  M  and  N  . 

The  allowable  virtual  processes  used  to  test  the  stability  of  the  system  can  be  visualized  as 
an  exchange  of  the  quantities  Xj^  between  the  subsystems  a  and  B  subject  to  the  constraints  of 
(2.5).   Thus,  a  virtual  process,  or  a  displacement  from  the  equilibrium  state,  has  the 
property  that 

6X°'  =   -  6X^   ,     i   =   1 r  (2.6) 

where  the  symbol  5  denotes  an  arbitrarily  small  displacement.   According  to  the  energy 
extremum  principle  (as  stated  by  Tisza),  any  virtual  process  which  takes  the  composite  system 
from  an  initially  unconstrained  equilibrium  to  a  constrained  equilibrium  must  result  in  a 
higher  energy  for  the  composite.   The  response  of  the  system's  energy  to  such  virtual 
processes  is  expressed  by  a  Taylor  series  expansion  about  the  initial  equilibrium  state,  i.e., 

=      5U"   .   i  6^   U"   ^   i   63    U"  +    .    .    .  (2.7) 

2  D 

.    5U^   .   i   6^    U^   .   1   6^    U^   .    .    .    . 
2  0 

where 

6U     =-        I       f-     6X.      -.        I        U.    Z.       ,  (2.8) 

1=1  1  1=1 

r  2  r 

6^U      =  I  ^^      6X      6X        =  I  U        Z   Z         .  (2.9) 

i,j  =  l       ^^^^j  '        ^  i,j  =  l         'J      '    ^ 

^^"    ^    .   \,     axixV    ^^i^^j^^k    ^    .   \^    ^jk^i^j^k    •  (2.10) 

i,J,k  =  l  1      J      k  -^  i,J,k  =  l 

"i   =  "3X7    '    '^ij    =    3X71X7    •    "ijk   =    3X.9X.9X,     '    ^""^   ^    =    ^^i    " 
1  1      J  1     J      k 

Substituting    (2.6,    2.8-10)    into    (2.7)    gives 


I      (u"  -  uh  z'^  .  i        I        (U".  ^  uf.)  Z°  Z" 

1=1  1112      ._j^^  ij  ij         i      J 


(2.11) 


a    S 
Since  a  and  B  were  initially  in  equilibrium,  U.  =  U.;  thus  the  first  order  form  is  identically 

zero.   The  quadratic  and  cubic  forms  are  simplified  by  the  following  consideration.   Since  the 
energy  is  represented  by  a  first  order  homogeneous  function, 

U"  [X^,  X^....,  X^]   =   k  U^  [X^^  X^ X^^]  (2.12) 

i^     ft 

where  k  is  the  ratio  of  the  sizes  of  two  systems  which  can  be  denoted  by  N  /N  . 
Also, 

X^  =  k  X^  .  (2.13) 

The  following  relationships  between  the  partial  derivatives  of  U  and  U   follow  from 
differentiating  (2.12)  and  making  use  of  (2.13): 


U"  =   U^   ,  (2.1H) 

1      1 


kU°'.   =   U^.   ,  (2.15) 


iJ     .  iJ 


k^  U";*.,   =   U^.,   .  (2.16) 


ijk      ijk 


Substituting  into  (2.11)  gives 


AU   ,     =  H^    I  U".  Z^Z'^  +  1^     I  U",  Z"Z"Z"  +  .  .  .   (2.17) 

system       2    ■    ";_.        13      ^  J  6       /-      ijk  1  j  k 

The  extremum  principle  requires  that  for  an  initially  stable  system,  all  virtual  processes 
must  result  in  the  change 

AU   ^    >   0   .  (2.18) 

system 

ot 
Hence,  the  quadratic  form,  which  dominates  the  expansion  as  the  virtual  displacements  Z.  are 

made  arbitrarily  small,  must  be  positive,  i.e., 

6^U  >   0   .  (2.19) 

In  (2.19)  we  divided  by  the  positive  constant  (1  +  k)/2,  and  dropped  the  superscript  a  since 
it  is  irrelevant  which  of  the  subsystems  is  considered. 

The  limit  of  stability  is  determined  by  the  condition 

6^U  =  0   .  (2.20) 


This  condition  is  of  particular  interest  in  this  discussion  because  a  critical  point  is  a 
stable  point  on  the  stability  limit.   The  stability  of  a  state  for  which  5^11  =  0  is  determined 
by  the  higher  order  forms.   Consider  the  cubic  form.   If  a  virtual  process  exists  such  that 
qSu  >  0,  then  because  the  cubic  form  is  an  odd  function  of  Z^,    there  also  exists  a  process 
(-Z^)  for  which  63u  <  0.   Therefore,  the  stability  of  a  critical  point  requires  that 

o^U  =  0     and     6^U  >   0   .  (2.21) 

In  a  pathological  case  where  6^U  =  0,  a  similar  argument  to  that  used  for  the  cubic  form  would 
require  that  a^U  =  0  and  6^U  >  0. 

In  fact,  it  is  normal  procedure  when  computing  critical  points  to  consider  only  the 
conditions  &^{i   =  0  and  &^U  =   0,  and  to  rely  on  experimental  evidence  to  infer  that  a  critical 
point  actually  exists  (i.e.,  is  stable). 

Before  moving  on  to  alternative  f orm.ulations  of  the  stability  conditioi'is,  there  is  one 
more  point  to  be  made  concerning  our  composite  system  model.   In  the  preceding  derivation  we 
found  the  relative  sizes  of  the  subsystems  a  and  S,  expressed  by  the  variable  k,  to  be 
irrelevant  to  the  stability  conditions.   There  is  a  special  value  of  k,  however,  which  adds 
valuable  intuition  about  the  kinds  of  virtual  processes  which  may  be  used  to  determine 
stability. 

Consider  the  case  where  one  of  the  subsystems,  B,  say,  becomes  very  large.   In  the  limit 
as  N^  ^  ",  k  ^  0,  and  by  eqs  (2.15-16)  the  second  and  higher  order  derivatives  of  U*^  go  to 
zero.   Recognizing  that  the  second  order  derivatives  of  U  are  the  same  as  the  first  order 

g 
derivatives  of  the  intensities,  U..  =  0  implies  that  the  intensities  of  6  remain  constant 

during  any  virtual  process.   This  corresponds  to  the  concept  of  ideal  reservoirs.   For 
example,  the  temperature  of  a  thermal  reservoir  remains  constant  when  heat  (entropy)  is  added 
to  or  removed  from  it.   Thus,  a  composite  system  consisting  of  a  subsystem  under  study  and  an 
ideal  reservoir  leads  naturally  to  the  consideration  of  virtual  processes  in  which  the 
intensities  are  held  constant. 

3.   Alternative  Formulations  of  the  Conditions  for  a  Critical  Point 

We  will  consider  a  critical  point  to  be  determined  by  the  conditions 

6^U  =  6^U  =  0   .  (3.1) 

Reformulations  will  be  considered  first  for  the  quadratic  form  and  subsequently  for  the  cubic 

form. 

3.1   The  Quadratic  Form 

A  stable  system  is  given  by  the  condition  6^11  >  0;  in  other  words  the  quadratic  form  is 
positive  definite.  We  shall  eventually  find  it  useful  to  express  the  results  concerning  the 
positive  definiteness  of  the  quadratic  form  in  terms  of  matrix  notations.   Hence,  the  notation 

is  introduced  here, 

6^U  =      z"^  U  Z   ,  (3.2) 


where  Z_   is  the  column  vector  (5S,  6V,  6ni  ,  ...,  6nj,_i)^, 


3^U 

U  is  a  symmetric  matrix  of  elements  U. .   =         , 

=  1 J  d  A  .  0  A  . 

1  J 

and  T  indicates  the  transpose  operation. 

A  common  technique  used  to  derive  conditions  which  determine  the  definiteness  of  a  quadratic 
form  is  to  diagonal ize  U  using  orthogonal  transformations.   However,  orthogonal 
transformations  do  not  turn  out  to  be  very  useful  in  thermodynamics  because  the  variables  do 
not  readily  admit  the  definition  of  a  metric  (see,  for  example,  V/einhold  [15]).   The  most 
useful  procedure  for  diagonalizing  U,  from  the  viewpoint  of  thermodynamics,  makes  use  of  the 
technique  completing-the-square.   The  diagonalization  of  U  and  subsequent  reformulations  of 
the  stability  conditions  will  now  be  considered  in  some  detail  for  a  binary  system. 

The  first  step  in  writing  the  quadratic  form  for  a  binary  system  as  a  sum  of  squares  is 

2  ^ 

-2,.     „    ,-       1  f„      „  „   ^  M      1 


where 


Continuing  to  complete  the  square, 


,(1)   -.2   r       .,,(1)  ^ 


6^U     - 


or 


where 


(3.3) 
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6^U   =   U^^  n^  *  U^2^  ^l   ^   U^^^  Z^  ,  (3.10) 
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^33      33    ,,(1)  ^^   ■* 
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Equation  (3.10)  implies  that  for  6^U   to  be  positive  definite,  the  conditions 

ui°^   =  "ii  >  0  .  (3.13) 

U^2^   >   0  ,  (3.n) 

U^^^   >   0  ,  (3.15) 

must  be  satisfied. 

The  first  step  toward  reformulating  the  stability  conditions  is  to  recognize  that  in  a 

stable  system  moving  toward  an  unstable  state,  the  first  quantity  which  will  go  to  zero,  and 

(2) 
thus  indicate  the  limit  of  stability,  is  U__  .   It  is  common  to  make  use  of  this  single 

condition  when  searching  for  the  stability  limit  of  a  system  of  interest  and  to  use  intuition 

gained  from  experimental  evidence  to  guarantee  that  the  search  is  started  from  a  stable  state. 

This  is  a  practical  approach.   The  purpose  of  indicating  the  implicit  assumptions  normally 

made  is  that  when  working  with  equations  of  state  and  Newton-Raphson  search  techniques  on  the 

computer,  we  can  unwittingly  project  the  solution  into  a  region  where  the  lower  order 

stability  conditions  are  not  met  and  hence  produce  incorrect  results. 

(2) 
To  prove  that  the  condition  U    is  the  first  to  go  to  zero,  consider  (3.6)  and  (3-12) 

which  indicate  a  general  recursion  formula  relating  successive  stability  conditions  and  can  be 
written  as: 

r  (k-2)  ^^ 
,,(k-l)     „(k-2)    ^  k(k-1)^  ,,  .,. 

\k     =  "kk    "7i^^2) •  (3.16) 

(k-l){k-1) 

(k-2)        (k-1  ) 
Under  the  assumption  that  the  system  is  initially  stable,  U,   >,,_,,  U     ,  and,  hence, 


U 


{k-2)  (k-2)  (k-2) 

..    are  all  positive.   If  we  let  'Jfi.-i  wu-i )  2°  to  zero,  then  assuming  that  U      remains 

,(k-2)  ,2  ^      ^    ^        ,,('<"'') 
J  ,    ,)   does  not  go  to  zero,  U 

K  ^  K   I  /  KK 

(k-2)  (k-1 ) 

have  become  zero  before  U,   ,,    ,  reaches  zero.   Therefore,  U      will  become  zero 


(k-2) 
before  U,   ,,    ,  does  as  the  stability  limit  is  approached. 

\K      i  J  \K      \   J 

If  we  assume  that  U    >  0,  and  use  the  recursion  formula  (3-16),  then  by  induction  the 

condition  which  will  go  to  zero  first,  and  thus  indicate  the  limit  of  stability,  will 

(r-1  ) 
be  U     ,  where  r  equals  the  number  of  rows  and  columns  in  U. 
rr 


(r-1  ) 
The  fact  that  the  stability  limit  is  determined  by  the  condition  U      =0  allows  us  to 

rr 

deduce  something  about  the  types  of  virtual  processes  which  may  be  used  to  test  stability.   In 
our  model  of  the  composite  system  we  have  imagined  virtual  processes  (Z)  as  the  exchange  of 
extensive  quantities  X^   between  the  subsystems  a  and  B.   The  diagonalization  of  the  quadratic 
form  for  the  binary  system  resulted  in  a  transformation  from  the  virtual  processes  (Z-],  Z2, 
Z3)  to  the  processes  ( m ,  ti2>  Z3).   In  general,  diagonalization  will  transform  the  processes 
(Zi,  Z2,  ...,  Zp)  to  the  processes  (ni,  112.  •••>  '^r-l  >  Zj,),   Since  the  limit  of  stability  is 
determined  by  the  coefficient  of  Zp,  the  virtual  processes  used  to  test  for  stability  need 
only  include  the  subset  given  by  (ni  =  0  (i  =  1,  ...,  r-1),  Z^) .      That  is,  we  may  arbitrarily 
set  the  variations  represented  by  the  m  to  zero.   Let  us  clarify  exactly  what  types  of 
processes  the  m  correspond  to. 

Consider  the  connection  between  variations  in  the  extensive  and  intensive  variables  for  a 
binary  system,  i.e., 
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(3.17) 


A  process  at  constant  temperature  implies  that 
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5T   =   0  , 


(3.18) 


or,  converting  notations 


"11  ^1  '   ^2  h  '   ^3  h     -     ° 


(3.19) 


and,  since  U-|  1  >  0  for  a  stable  system,  we  divide  by  it  to  obtain 
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0  . 


(3.20; 


Comparing  (3.20)  and  (3-5)  leads  us  to  the  conclusion  that  m  =  0  is  equivalent  to  a  virtual 
process  at  constant  temperature. 

Since  the  stability  limit  may  be  determined  by  using  constant  temperature  processes,  it 
is  natural  to  reformulate  the  quadratic  form  in  terms  of  the  Helmholtz  free  energy,  i.e.. 


6h 


I  A. ,  Z.Z. 


(3.21) 


where 


Diagonalizing  leads  to 


(6V,  dn^)"^  .  (3.22) 


2  2     (1)2 

6^A   =   A^^  n^  +  A^2   ^2  .  (3.23) 

And  the  limit  of  stability  is  determined  by 

(1  )  12 

4    =   ^22  -17;-   =   °  •  (3.2^) 

The  stability  limit  takes  a  simpler  form  in  terms  of  A  because  the  constant  temperature 
constraint  is  implicit  in  A  which  is  the  first  Legendre  transform  of  U  with  respect  to  S. 

Let  us  carry  the  present  line  of  thought  one  step  further.  Consider  virtual  processes 
which  occur  at  constant  temperature  and  pressure.   The  constraints  placed  on  the  variations  of 
the  extensive  variables  are 

^1  h    '   "12^2^  U^3Z3   =   6T  =  0  , 
^21  h    ^   "22  ^2  ^  "23  ^3   =  -«P   =  0  • 


Eliminating  Z-|  from  these  equations  gives 

Z   + y Z   =  0  .  (3.25) 

"12 

U   -  —-=- 

22    U^^ 

Comparing  (3.25)  and  (3.11)  shows  that  n2  =  0  is  equivalent  to  a  process  at  constant 
temperature  and  pressure.   Hence,  the  stability  conditions  can  be  expressed  in  the  Gibbs  free 
energy  formulation,  i.e., 

6^0  =   G^^  Z^  ,  (3.26) 

where 

^1   =  ^"1  • 

The  limit  of  stability  is  given  by  the  condition 

G^^   =  0  .  (3.27) 

For  multicomponent  mixtures  we  can  generalize  the  conditions  for  the  limit  of  stability 
as  follows.   In  the  energy  representation, 

U^^'^^   =  0  .  (3.28) 

rr 


Considering  virtual  processes  at  constant  T  leads  to 

(r-l)(r-1)     "  • 


(3.29) 


and  considering  processes  at  constant  T  and  p  leads  to 

^(r-2)(r-2)   "   ° 


(3.30) 


Beyond  this  point  we  run  out  of  commonly  used  symbols  to  represent  the  Legendre  transforms  of 
U,  hence,  the  notation  [U^]  is  introduced  to  represent  the  i^^   Legendre  transform  of  U. 
Considering  processes  at  constant  T,  p,  p-|  leads  to 

.  (r-4) 
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(3.31) 


(r-3)(r-3) 


and  considering  processes  where  r-1  intensities  are  held  constant,  i.e.,  constant  T,  p,  \i-\  , 
.  . . ,  Vq-2 I  leads  to 

(0) 


r-1 


[U^"'] 


0 


(3.32) 
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There  are  r  equivalent  formulations  of  the  stability  conditions  for  the  particular 
ordering  of  independent  variables  we  have  chosen.  Many  more  formulations  are  possible  by 
taking  different  permutations  of  the  order  S,  V,  n-]  ,  ...,  nQ_i  .   Of  all  possible  formulations, 
there  is  really  only  one  which  interests  us  here;  that  is  the  formulation  in  terms  of  the 
Helmholtz  free  energy  because  most  equations  of  state  are  written  with  T,  V,  n-|  ,  ...,  n^  as 
the  independent  variables. 

In  addition  to  reformulating  the  stability  conditions  in  terms  of  the  Legendre  transforms 
of  the  energy,  there  is  another  important  reformulation  which  we  consider  now.  As  before,  the 
motivation  will  come  from  working  with  the  specific  example  of  the  binary  system. 

Consider  the  stability  conditions  given  in  eqs  (3.13-15).   The  first  condition  can  be 
trivially  rewritten 
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(0) 

11 
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11 
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:3.33) 


where  A-|  is  the  first  discriminant  of  the  matrix  U.   The  m^'^  discriminant  of  a  matrix 

is  the  determinant  of  the  submatrix  obtained  by  deleting  all  elements  which  do  not 

simultaneously  lie  in  the  first  m  rows  and  columns  (see,  for  example,  Hildebrand,  p.  52). 
The  second  condition  can  be  rewritten. 


U^^  ^   =  — ^  (U   U   -  U^  ) 
22      U    ^11   22    12-^ 


(3.34) 


or,  due  to  the  symmetry  of  U, 
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The  third  condition  can  be  rewritten. 
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The  quadratic  form  for  a  binary  system  can  then  be  written, 


and  the  condition  for  the  stability  limit  translates  into 


For  the  binary  system,  U  is  a  3  x  3  matrix,  hence  (3.40)  becomes 


In  general. 


(3.36) 


(3.37) 


(2)     ^3 
"33   =  A^-  (3.38) 


2  2    "^2   2    ^3   2 

5  U  =  A^  n^  .-  n^  -^  Z3  ,  (3.39) 


A^  =  0  .  (3.40) 


A   =  det  j  U  I   =   0  .  (3.41) 


1       2  r-2        r-1 


and  the  stability  limit  can  be  represented  by 

U^^'^^   =   A^  =   det  I  y  I  =  0  .  (3.43) 

The  condition  for  the  stability  limit  expressed  as  Ap  =  0  is  in  the  final  form  as  given 
by  Gibbs  and  is  in  the  form  used  for  computation  by  Heidemann  and  Khalil.   Equivalent  forms 
are  obtained  when  working  with  a  Legendre  transformation  of  the  energy  representation.   In 
general,  the  stability  limit  can  be  expressed 

.  (r-i-1) 
[U^]  =   A^  .   =   0  .  (3.44) 

(r-i)(r-i) 

Before  moving  on  to  a  discussion  of  the  cubic  form,  there  is  one  final  point  to  be  made 
regarding  the  quadratic  form  which  is  suggested  by  the  condition  for  the  limit  of  stability 
written  as  a  determinant.   Consider  the  energy  formulation.   At  the  stability  limit, 
eq.  (3.43)  implies  that  U  is  a  singular  matrix.   We  have  mentioned  previously  that  U 
represents  a  transformation  between  the  extensive  and  intensive  system  variables.   When  a  sys- 
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tern  is  in  a  regular  state  (i.e.,  U  is  not  singular)  then  a  set  of  r  extensive  variables 
uniquely  specifies  a  set  of  r  intensive  variables  (recall  that  x^+i  has  been  fixed  to  specify 
the  size  of  the  system).   This  is  in  accord  with  Gibbs'  phase  rule,  i.e.,  the  number  of 
independent  intensive  variables  in  a  single  phase  system  is  e+1  =  r.   At  the  limit  of 
stability,  the  rank  of  U  is  reduced  from  r  to  r-1 ;  hence,  there  are  only  r-l  independent  in- 
tensive variables.   Since  U  is  singular,  there  exists  a  non-trivial  solution  to  the  homogene- 
ous equation 

y  Z  =  0  .  (S.'JS) 

Thus,  the  virtual  process  which  causes  the  quadratic  form  to  be  zero  at  the  stability  limit 
defines  a  unique  direction  in  Gibbs  space  within  an  arbitrary  sign;  this  process  is  referred 
to  as  a  critical  displacement.   Furthermore,  the  critical  displacement  expressed  by  eq  (3.^5) 
represents  a  virtual  process  in  which  the  r-1  independent  intensities  are  held  constant 
(recall  eq  (3.17)). 

3.2  The  Cubic  Form 

The  condition  6->U   =  0  is  necessary,  but  not  sufficient,  to  prove  the  stability  of  a 
system  on  the  stability  limit.   In  this  section  two  formulations  of  the  condition  6^U  =   0  will 
be  presented.   Analogous  formulations  in  terms  of  the  Legendre  transforms  simply  involve  fewer 
independent  variables. 

The  previous  section  established  that,  for  a  system  on  the  stability  limit,  there  exists 
a  virtual  process,  unique  within  an  arbitrary  multiplicative  constant,  called  the  critical 
displacement  for  which  the  quadratic  form  becomes  identically  zero.   Therefore,  the  critical 
displacement  is  the  only  process  which  needs  to  be  considered  when  evaluating  the  cubic  form. 

The  cubic  form  is  a  special  case  of  a  general  tri-linear  form.   Its  linearity  allows  it 
(2.10)  to  be  written  as 

^'"     -      X      h-i-     t,    ^      WIT     ^i^J^  ^3.46) 

k=l  k        i,j=1  1      J 

or 

&h     =       I       \^     CS^U]    .  (3.^*7) 

k=1  k 


or 


6^U     =      Z-V   (6^U)    ,  (3.48) 
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where  V  is  a  gradient  operator  with  components  -^r—  .   The  operator  Z»V  produces  a  directional 


derivative,  that  is,  the  derivative  in  the  direction  of  the  vector  1_   in  Gibbs  space.   (Gibbs 
space  is  taken  here  to  be  spanned  by  the  variables  S,  V,  n-|  ,  ...,  n^--^.)     Hence,  eq   (S.'iS) 
indicates  that  the  cubic  form  necessary  to  show  the  stability  of  a  system  on  the  stability 
limit  is  the  derivative  of  the  quadratic  form  in  the  direction  of  the  critical  displacement. 
To  evaluate  the  cubic  form,  the  direction  of  the  critical  displacement  vector  must  be  known; 
this  vector  is  determined  for  systems  on  the  stability  limit  by  solving  eq  (3-^5). 

Gibbs  formulated  the  second  condition  for  a  critical  point  in  a  way  which  avoided  having 
to  explicitly  determine  the  critical  displacement  vector.   His  formulation  is  more  general 
than  is  the  specification  that  S^U  =  0,  but  it  contains  a  redundancy  which  is  extremely  costly 
from  a  computational  point  of  view.   To  formulate  Gibbs'  second  condition  we  proceed  as 
follows. 

The  critical  displacement  is  one  in  which  r-1  of  the  r  intensities  are  held  constant.   To 
consider  virtual  processes  of  this  type,  choose  any  r-1  of  the  r  equations  represented  by 
eq  {3-^5)   to  act  as  constraints  on  the  r  independent  extensive  variations  1^,    i  =  1,  ...,  r. 
For  definiteness  choose  the  first  r-1  equations.   These  equations  impose  the  constraint  that 
3U 


the  r-1  intensities  -i-rr-    ,    i 

oX  . 

1 


1 r-1,  are  constant  in  any  virtual  process.   The  only 

independent  variation  is  6Xp  which  may  be  chosen  to  be  arbitrarily  small  as  required  to 
determine  local  stability.   For  the  final  equation  consider  an  analog  to  the  cubic  form;  that 
is,  rather  than  take  a  directional  derivative  of  the  quadratic  form,  differentiate  the 
discriminant  A^  which  is  an  equivalent  indicator  of  the  stability  limit.   Just  as  Z'V(62u) 
must  be  zero  for  a  critical  point,  so  must  Z_'V(Ap).   Then  our  final  set  of  equations  is 
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(3.iJ9) 


A  nontrivial  solution  to  eq   (3.^9)  requires  the  matrix  to  have  a  zero  determinant,  which  is 
Gibbs'  second  condition  for  a  critical  point. 
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3.3  Computational  Considerations 

Let  us  first  consider  the  computations  required  by  Gibbs'  formulation  for  a  critical 
point.   The  conditions  are  that  Ap  =  0  (eq  (S.'^S))  and  that  A^  =  0,  where  A^  is  the 
determinant  of  the  matrix  in  eq   (3- '^9).  Both  conditions  involve  r^^   order  determinants, 
however,  the  second  condition  contains  derivatives  of  A^  within  the  determinant  A^.   An 
exact  formula  for  the  derivative  of  a  determinant  is  derived  by  expressing  the  determinant  as 
a  multiple  matrix  U,  then  (Courant  and  John  [1]) 

A   =   det  I  U  I  =  U,  X  U^  X  . . .  X  U    -U  (3.50) 

r        1=1     -1-2        -r-1  -r 

where  x  is  a  vector  cross  product  and  •  is  a  vector  dot  product.   Differentiating  (3.50)  gives 

3A^     3U^ 

^  =   33^  X  U^  X  ...  X  U^_^.U^ 
1       1 

9U 
.   U^  X  3^  X  ...  X  U^_^.U^  (3.51) 

1 


3U 
.   U^  XU^X  ...  XU^_^.^  ; 

1 

hence,  the  derivative  of  an  r^^   order  determinant  equals  the  sum  of  r  determinants,  in  each  of 
which  one  column  is  replaced  by  its  derivative  with  respect  to  the  appropriate  independent 
variable.   Hence,  to  evaluate  Gibbs'  conditions  for  a  critical  point  of  a  single  state 
requires  the  evaluation  of  (r^  +  2)  r^^   order  determinants.   The  most  efficient  way  to  compute 
a  determinant  is  to  factor  the  matrix  and  then  compute  the  product  of  the  pivots.   Matrix 
factorization  requires  approximately  r3/3  operations;  therefore,  the  computational  effort 
required  to  evaluate  Gibbs'  conditions  is  proportional  to  r^. 

The. efficiency  of  the  method  proposed  by  Heidemann  and  Khalil  for  computing  critical 
points  comes  from  their  realization  that  if  eq   (3.^5)  is  solved  for  the  critical  displacement 
vector  Z,    then  the  cubic  form  may  be  evaluated  directly.   Direct  evaluation  of  the  cubic  form 
requires  r3  operations.   It  will  be  shown  later  that  for  some  equations  of  state,  such  as  the 
general  two-constant  cubic  equation  to  be  considered,  the  cubic  form  can  be  analytically 
summed  in  such  a  way  as  to  reduce  the  computational  effort  required  for  its  evaluation  to 
being  proportional  to  r^.  When  the  cubic  form  may  be  simplified  in  that  way,  the  major 
portion  of  the  critical  point  calculation  is  required  to  factor  the  matrix  U;  thus  the  overall 

effort  remains  proportional  to  r3. 

In  addition  to  the  difference  in  computational  effort  required  by  the  two  formulations 
presented  above,  they  require  distinct  iteration  schemes  which  will  be  discussed  here.  The 
conditions  of  criticality  will  be  expressed  in  terms  of  the  Helmholtz  free  energy,  and  the  set 

of  independent  variables  is  taken  as  (T,  V,  n-] n^-i).   Since  there  are  two  conditions  of 

criticality,  we  will  specify  the  composition  of  the  mixture  of  interest  and  then  solve  two 
equations  for  the  values  of  T,  V  at  a  critical  point.   Critical  pressure  is  subsequently 
obtained  from  the  equation  of  state. 
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Gibbs'  formulation  results  in  the  equations 

Vl   =   °   •      ^^-1   =   °   •  ^3.52) 

These  equations  may  be  solved  using  a  Newton-Raphson  iteration  scheme.   The  derivatives  of  the 
determinants  required  for  the  iteration  are  normally  computed  numerically.   Equation  (3-52) 
may  be  evaluated  at  any  state  point,  and  the  critical  point  results  will  be  correct  as  long  as 
the  stability  limit  is  approached  from  a  stable  region. 

Heidemann  and  Khalil  formulated  the  critical  point  as 

A^_^   =  0   ,     6^A  =  0   .  (3.53) 

The  solution  procedure  must  be  modified  from  a  straightfoward  Newton-Raphson  technique  because 
the  cubic  form  can  be  evaluated  only  at  a  state  on  the  stability  limit.   Therefore,  if  the 
volume  is  specified,  the  temperature  on  the  stability  limit  may  be  determined  by  solving 
A)---]  =  0  using  a  one-dimensional  Newton-Raphson  iteration.   Then  the  critical  displacement 
vector  is  determined  by  solving  A  Z_  =  0,  and  the  cubic  form  is  evaluated.   The  search  has  now 

been  reduced  to  a  one-dimensional  search  for  the  critical  point  along  the  limit  of  stability. 
The  derivative  of  the  cubic  form  with  respect  to  volume  along  the  stability  limit  is 
determined  numerically,  and  a  Newton-Raphson  Iteration  employed  to  update  the  volume.   This 
procedure  of  nested  one-dimensional  searches  is  continued  until  a  critical  point  is  found.   In 
addition  to  the  efficiency  of  computing  &^k   rather  than  A^_-|  ,  the  formulation  of  Heidemann 
and  Khalil  appears  to  have  a  larger  radius  of  convergence.   Heidemann  and  Khalil  reported  that 
a  single  strategy  for  the  initial  guess  was  sufficient  to  compute  critical  points  for  the 
large  group  of  mixtures  they  considered. 

iJ.   The  [lethod  of  Heidemann  and  Khalil 

Heidemann  and  Khalil  chose  to  work  with  the  stability  conditions  written  in  terms  of  the 
Helmholtz  free  energy;  this  is  the  logical  choice  when  working  with  equations  of  state  with  T, 

V,  n-i n^  as  the  independent  variables.   Thus,  all  virtual  processes  are  carried  out  at 

constant  temperature.   Furthermore,  Heidemann  and  Khalil  fixed  the  size  of  the  system  by 
arbitrarily  holding  V  constant,  rather  than  N  as  was  done  in  the  theory  sections  of  this  note; 
this  results  in  the  independent  variable  set  ni ,  n2,  ...,  n^,  where  c  is  the  number  of 
components  in  the  mixture.   Their  claim  that  this  choice  of  variables  "has  the  effect  of 
producing  symmetrical  quadratic  and  cubic  forms"  is  misleading.   The  symmetry  of  the  quadratic 
and  cubic  forms  is  predicated  on  the  order  of  differentiation  in  the  second  and  third  order 
derivatives  of  the  free  energy  being  arbitrary.   In  general,  the  variables  n-j  ,  n2.  ....  n^  are 
convenient  from  the  point  of  view  of  computations  because  all  the  derivatives  of  the  free 
energy  will  have  the  same  form.   Furthermore,  for  the  particular  equation  of  state  to  be 
examined  in  this  work,  this  choice  of  variables  results  in  significant  simplifications  in  the 
analytical  expressions  of  the  quadratic  and  cubic  forms. 

In  terms  of  the  Helmholtz  free  energy  (A),  the  stability  limit  is  expressed  by  the 
condition  (dropping  the  subscript  on  Ar>_i  ) 

A  s  det  I  Q  I  =  0  (1.1) 

where  the  components  of  the  c  x  c  matrix  Q  are 
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1  J 

For  a  specified  value  of  V,  eq  (4.1)  is  used  to  determine  the  value  of  T  on  the  stability 
limit.  Because  of  the  expense  of  computing  -r—  analytically,  and  because  the  computation  of  A 

o  1 

itself  is  the  most  time  consuming  part  of  the  critical  point  computation,  a  quasi-Newton 

iteration  is  used  to  solve  [H.]).      This  technique  involves  computing  -r;r   numerically  on  the 

oT 

first  iteration,  i.e., 

9A(T  )      A(T   +  6T)  -  A(T  ) 

^  =  °     °-  '  (1)  3) 

3T  6T  ^    '^' 

where  6T  =  10~°  Tq,  Tq  is  the  initial  temperature  guess,  and  A(T)  indicates  the  determinant 
evaluated  at  T  (V,  n-\  ,    .  , . ,  n^  are  fixed  during  the  temperature  iteration).   On  subsequent 
iterations  use 

3A(T.)      A(T.)  -  A(T.  , ) 

'-'-  (4. it) 


3T         T.  -  T.  , 
1    1-1 


th 


where  Tj^  is  the  temperature  after  the  i^''  iteration  which  is  obtained  from 

A(T.) 


3T 

The  quasi-Newton  iteration  exhibits  superlinear  convergence,  and  requires  only  one  evaluation 
of  A  per  iteration  (after  the  first).   The  determinant  computation  is  performed  by  subroutines 
DSPFA  and  DSPDI  from  the  LINPACK  software  (Dongarra  et  al .  [2]).   The  first  routine  factors  a 
symmetric  indefinite  matrix  (stored  in  packed  form)  and  the  second  computes  the  determinant  as 
a  product  of  the  pivots. 

Once  the  iteration  to  find  the  stability  limit  has  converged,  the  cubic  form  is 
evaluated:   the  first  step  in  this  process  is  to  determine  the  critical  displacement  vector  Z; 
i.e.,  solve 

Q  Z  =  0   .  (4.6) 

In  general  this  step  is  efficiently  performed  by  the  inverse  iteration  technique  (Peters  and 
Wilkinson  [11]).   In  the  implementation  described  here,  Z  is  obtained  from  the  LINPACK 
subroutine  DSPCO. 

Because  Z   is  determined  only  within  an  arbitrary  multiplicative  constant,  it  must  be 
subjected  to  further  constraints  to  avoid  arbitrary  changes  in  the  sign  or  magnitude  of  the 
cubic  form  evaluated  along  the  stability  limit;  such  changes  falsely  indicate  the  location  of 
critical  points.   The  "length"  of  Z  is  normalized  by  the  constraint 

z"^  Z  =  1   .  (4.7) 
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We  determine  the  direction  of  !_   by  specifying  that  its  first  component  Snj  be  positive  at  the 
initial  guess  of  the  volume.  Constraining  (>n^   to  be  positive  to  fix  the  direction  of  !_   cannot 
be  continued  as  we  move  along  the  stability  limit  searching  for  the  zero  of  the  cubic  form. 
Because  it  is  possible  for  6ni  to  go  through  a  zero  as  Z_   changes  with  V  along  the  stability 
limit,  constraining  6ni  to  always  be  positive  would  be  equivalent  to  changing  the  direction  of 
Z  at  a  value  of  the  volume  where  6ni  goes  through  a  zero.   The  technique  used  to  ensure  that 
the  direction  of  Z_  varies  continuously  as  a  function  of  V  along  the  stability  limit  is  the 
following.  The  null  vector  from  the  previous  iteration  on  V  is  saved  and  projected  onto  the 
null  vector  computed  for  the  current  iteration.   If  the  projection  is  positive  we  assume  that 
the  direction  of  Z_  is  changing  continuously.   This  final  condition  on  Z_   may  be  written 

z'[_^  Z .  >  0  ,  (4.8) 

where  Zj^  is  the  null  vector  from  the  i^'^  iteration  on  V.   If  Z_]_   does  not  satisfy  C^.S),  it  is 
multiplied  by  -1  so  that  it  does. 

Michelsen  and  Heidemann  [7]  have  recommended  that  discontinuities  in  Z  be  avoided  by 
making  use  of  the  condition 

z"^  Z.  >  0   ,  (1.9) 

where  Iq   corresponds  to  the  initial  guess  Vq.   Depending  on  how  far  Vq  is  from  the  critical 
point,  it  is  quite  possible  that  Z  could  change  enough  along  the  stability  limit  that 
eq  (M.9)  would  result  in  an  arbitrary  change  in  direction.   It  seems  safer,  and  no  more 
costly,  to  employ  eq  (4.8). 

The  summation  of  the  cubic  form  can  in  general  be  reduced  from  c3  terms  to  c(c+l  ) (c+2)/6 
by  making  use  of  its  symmetry  properties,  i.e., 

c  c   c   c 

6'A  =     I  A..,  Z.Z.Z,  =   y    y    y  h..,  a..,  Z.Z.Z,  (4.10) 

.  ■  ,  ,   ijk   1  J  k   ....  .  ,   ijk   ijk   1  J  k 
i,J,k  =  l    ^     ^     k=j  j  =  i  1  =  1    -^    -^     ^ 

where 

1;  i=j=k 
hijk  =     3;  i=j^=k,  or  i  =  k|=j,  or  i=j=j=k 
6;  i=t=j=fk  and  i=|=k 
For  the  special  case  of  the  two-constant  cubic  equation  of  state  to  be  considered  here,  the 
cubic  form  has  been  simplified  so  that  only  a  double  summation  is  required  for  its  evaluation. 
Michelsen  [6]  suggested  that  the  cubic  form  could  be  evaluated  numerically  because  it  is 
the  derivative  of  the  quadratic  form  in  the  direction  of  the  critical  displacement  vector; 
hence, 

-        z'^Q(T,V,n  +  £Z)Z  -  z'^Q(T,V,n)Z 
6^A  S  C  =  TT— ZTT (4.11) 


¥ 


where         n  is  a  vector  of  the  mole  numbers,  i.e.,  (n-],  n2 Hq)', 

e  is  a  small  number,  i.e.,  0(10~6)^ 
zT  Q(T,V,n)Z  =  0  on  stability  limit, 

lleZjl    =    e||z)|    =   £,    since  ||z||    =   Z^Z   =    1. 
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This  expression  for  the  cubic  form  is  particularly  useful  when  the  required  derivatives  of  A 
are  not  expressible  analytically. 

Once  the  stability  limit  is  located  and  the  cubic  form  is  evaluated  at  the  current  value, 
Vj^,  the  volume  is  updated  by 

C(V.) 
^i.1  =^i  -C^V  '  ^'^•^2) 

where  C'(Vj^)  represents  the  derivative  of  C  with  respect  to  the  volume  (dC/dV)  or  some 
estimate  of  it.   The  original  method  of  Heidemann  and  Khalil  approximated  dC/dV  numerically, 
i.e., 

d  C(V.)    C  (V.  +  6V)  -  C  (V.) 

dV  6V  ^^'':>J 

where  6V  =  lO'^Vj^.   This  approximation  is  expensive  because,  when  the  volume  is  incremented, 
eqs   (^.1)  and  (M.6)  must  be  solved  again  before  evaluating  C{y^   +   6V).   A  good  approximation 
can  be  made  for  T(Vj^  +  6V)  with  only  an  inexpensive  calculation  (see  below)  so  that  only  one 
or  two  iterations  are  required  to  solve  (4.1).   But  we  must  keep  in  mind  that  the  single  most 
expensive  calculation  in  this  method  is  the  factorization  of  Q;  thus  we  should  try 
to  minimize  the  number  of  evaluations  of  the  cubic  form.  Michelsen  and  Heidemann  [7] 
presented  a  quasi-Newton  method  which  requires  only  a  single  evaluation  of  C  for  each 
iteration  on  V.   In  general,  they  used 

C(V.)  -  C(V.   ) 

C'(V.)  =  y  \  y ^-^  .  (4.14) 

i    i-1 

Two  values  of  C  are  needed  to  start  this  iteration.   One  option  Is  to  use  (4.13)  for  the  first 
iteration,  then  switch  to  (4.14)  for  subsequent  iterations.   Michelsen  and  Heidemann  chose  to 
use  a  modified  form  of  (4.13)  for  the  first  iteration.   To  evaluate  C(Vj^+6V)  they  make  two 
approximations:   first,  that  Z^i  =  Zq,  and  second,  that  T  (Vj^  +  6V)  is  obtained  from  a  Taylor 
series  expansion,  i.e., 

T(V.  +  6V)  =  T(V.)  +  SV  IJ  I     ,  (4.15) 

i 

where  T(Vj^+6V)  is  the  first  order  approximation  to  the  temperature  on  the  stability  limit  at 
(Vj^+6V),  T(Vj^)  is  the  temperature  on  the  stability  limit  at  Vj^,  and  dT/dV  is  the  derivative 
along  the  stability  limit.   To  obtain  an  expression  for  this  derivative,  consider  the 
stability  limit  for  a  given  mixture  (fix  n)  to  be  parameterized  in  terms  of  V.  Since  the 
stability  limit  is  defined  by  Z^Q  Z^  =  0,  then 

ri    T       '^^  T  '^S      T  '^- 

dv^i5^^=dV   5Z^Z   d^^^^Q-=0 

(4.16) 
T  dQ 

=  ^  d?^=° 
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Since  Q  Z  =  Z^Q  =  0.   By  chain  rule, 


dQ   3Q   3Q  .„ 

dV   9V   3T  dV  ^■^•'  1 1 


Substitution  of  (^.17)  into  (M.16)  gives 

dV  -      3Q   •  ^''•^^^ 

^  3f^ 

These  two  approximations  allow  an  estimate  of  C(Vj^+6V)  which  requires  no  additional 
factorizations  of  Q.  After  the  initial  iteration  on  V,  C(Vi)  in  (M.14)  is  evaluated  exactly, 

and  the  initial  guess  of  temperature  at  a  new  volume  V^  is  made  analogously  to  (4.15),  i.e., 

dT 


T  (V.)  =  t[v.  J  +  (v.  -  V.  J  -Si 
0  1     ^  1- V    ^1    1-1   dV 


(i|.19) 
^i-1 


Although  the  quasi-Newton  method  requires  only  half  the  number  of  evaluations  of  the 
cubic  form  per  iteration  on  V  as  does  the  full  Newton-Raphson  technique  described  initially, 
its  convergence  is  slower  (superlinear  rather  than  quadratic).   In  general,  the  cost  of 
computing  critical  points  by  various  methods  will  depend  on  the  cost  of  computing  the 
derivative  C'(Vj^)  relative  to  the  convergence  rate  which  can  be  attained  with  that 
approximation.  This  must  be  determined  by  numerical  experiments. 

4.1   Analytical  Expression  for  the  Volume  Derivative  of  the  Cubic  Form 

It  may  be  advantageous  to  retain  the  quadratic  convergence  of  a  Newton-Raphson  iteration 
if  the  derivative  of  the  cubic  form  with  respect  to  volume  can  be  computed  more  efficiently. 
In  this  section  is  presented  an  analytical  derivation  of  dC/dV  which  requires  one  additional 
factorization  of  an  augmented  Q  matrix  to  be  computed  for  each  volume  iteration.  The 

procedure  is  substantially  more  efficient  than  the  forward  difference  approximation  of  dC/dV 
(eq  4.13)  which  requires  two  evaluations  of  C  at  each  volume  iteration. 

We  begin,  as  previously  for  dT/dV,  by  considering  the  stability  limit  to  be  parameterized 
in  terms  of  V.   Then 

c      dA .  dZ . 

^=   I  [— ^  Z-Z-Z,  +3A..  -;^  Z.Z,  ]  (4.20) 

dV   .  .  ,  ,     dV    1  J  k     ijk  dV   j  k-' 
i,j,k=l 

where 


dA .  . ,     3A .  . , 
ijk  _   ijk 

dV 


(4.21) 


The  partials  OA^j^-ZSV)  and  (3Aij|</3T)  may  be  obtained  by  differentiation  of  the  expression 

dZ 
for  Aj^ji^  obtained  from  the  equation  of  state.  The  derivative  —  is  obtained  from  the 

following  considerations.   Differentiating  the  expression  Q  Z,  which  is  zero  on  the  stability 

limit,  leads  to 
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dZ 
dV 


dQ 
d^^ 


(4.22) 


dZ 
This  expression  could  be  solved  for  —  except  that  Q  is  singular  on  the  stability  limit. 

dZ 

Therefore  we  must  find  another  condition  on  ^rr  which  will  be  used  to  replace  the  row  in  Q 

dV  = 

containing  the  zero  pivot.  The  appropriate  condition  is  from  a  simple  result  of  vector 
calculus:   the  derivative  of  a  vector  of  constant  length  is  orthogonal  to  that  vector,  i.e., 


-  dV 


(4.23) 


This  condition  can  be  combined  with  (4.22)  and  the  result  written  as  the  solution  of  c+1 
equations  in  c  unknowns, 


' 

" 

dQ 

Q 

dZ 

dV^ 

vT' 

dV~ 

= 

______ 

_ 

L 

(4.24) 


5.   Implementation  for  a  Two-Constant  Cubic  Equation  of  State 

Expressions  for  the  derivatives  of  the  Helmholtz  free  energy  are  required  to  evaluate  A, 
C,  and  their  derivatives  with  respect  to  T  and  V.   These  expressions  are  given  here  for  a 
general  two-constant  cubic  equation  of  state  (TCC  EOS)  following,  for  the  most  part,  the 
notation  of  Michelson  and  Heidemann  [7]. 

The  TCC  EOS  is  written 


NRT 
V-Nb 


N^a 


(V  +  6  Nb)(V  +  &^Hb)    ' 


(5.1) 


where 


—rr        )        n.n  .a.  . 

n2  i.L  '  J  ^j 


(5.2) 


1=1 


(5.3) 


N  =   I  n.   , 
1=1 


a..  =  ^..  (a.,  a..) 


1/2 


(5.4) 


a. .  =  n  R^T^./P  .  fl  +  m.[l 
11    a    ci  ci  '     1 


/2i  i2 


(T/T  .)''^] 
ci 


(5.5) 


m .  =  m(a).  )   , 
1      1 
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b .  =  n^  RT  . /P  . 
1    b   ci  ci 


(5.6) 


For  the  SRK  equation  of  state, 


6,  =  1,  6„  =  0,  Q  =  0.M2748,  n  =  0.0866iJ   , 


(5.7) 


m(u.)  =  0.480  +  1.574u.  -  0.176a)/ 


(5.8) 


and  for  the  PR  equation  of  state, 


6,  =  1  +  /2,  6,  =  1  -  /2,  n  =  0,45724,  Q     =  0.07780   ,  (5.9) 

I  ^  cL  D 


m(a).)  =  0.37464  +  1.54226  w.  -  0.26992  w.' 


(5.10) 


Heidemann  and  Khalil  chose  to  proceed  from  the  expression  for  fugacity 


RT  In  f.=  - 

1 


-3P 

■Sn. 

1 


RTn  "i  ^^ 

^]  dV  +  RT  in  ^;j— 


T,V,n 


(5.11) 


Their  choice  of  (n-|  ,  n2.  ...,  nj,)  as  the  independent  variables  results  in  all  derivatives 
containing  differentiations  with  respect  to  a  mole  number.   However,  if  one  of  the  mole 
numbers  in  the  independent  variable  set  were  to  be  replaced  with  the  volume,  then  one  element 

of  g  (i.e.,  — -)   would  contain  no  differentiations  with  respect  to  mole  number.   In  that  case, 

3V 
those  derivatives  would  have  to  obtained  by  differentiating  the  Helmholtz  free  energy  rather 

than  the  fugacity.  To  be  consistent  with  any  derivatives  obtained  from  the  fugacity,  the  free 

energy  must  be  relative  to  a  reference  state  of  pure  ideal  gases  at  the  system  temperature  and 

unit  pressure,  i.e.  (see  Prausnitz  [12],  p.  40), 


A  -  A 


CO  1=1  1 


(5.12) 


where 


^°=  I 


i  =  1 


A.  =  Helmholtz  free  energy  of  pure  i  in  an  ideal  gas  state  at  the  system  T,  and  p  =  1 


When  this  expression  is  differentiated  with  respect  to  mole  numbers, 


3  A. 
: 

3  n. 


y?  -  RT 


(5.13) 


T,V 
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where  y  =  the  chemical  potential  of  pure  i  in  an  ideal  gas  state  at  the  system 
i 

T,  and  p  =  1 . 

The  Helmholtz  free  energy  for  the  TCC  EOS  is 

c        n.RTF 
A-A°  =  RT  l^    n.    in^j^-^F^  ,  (5.14) 

where  Fj^  are  functions  of  the  dimensionless  volume  K  {=   V/Nb)  which  are  summarized  at  the  end 
of  this  section.   Differentiating  (5.14)  with  respect  to  mole  numbers  gives 

n.RTF  e. 

RT  in  f.  =  RT  Un  -\^  ^    B.  F^  ]  -  |  (,.  p^  .  _i  F^]  ,  (5.15) 


where 


c  n .  a.  . 

j=i 

b. 

The  elements  of  Q  are  obtained  from  further  differentiation  of  (5.15)  with  respect  to  mole 
numbers ,  i.e., 

9  In  f .  N  6.  . 

f^^ij  =  ^^^^-37^  — ^T,V,  n,  =  ^^t^^  ^  ^^i  ^  2j)  ^1  ^  ^i^j  ^1  J 
J         k         1 


(5.18) 


where  6j_j    is  the  Ki^onecker  delta,  and  n^.  held  constant  in  the  differentiation  signifies  that 
all  mole  numbers  are  held  constant  except  nj . 

To  evaluate  the  cubic  form  one  further  differentiation  of  (5.18)  with  respect  to  the  mole 
numbers  is  required,  i.e., 

?       ^     3^  in  f 
'   ^Jk  =  "   «^  (3?r-3-n-^T,V,  n^ 

=  RT[-  ii-iii  .  (6.B.  .  S.B^  ^  6.6,)F^2  ^  ^  e.B.B,  F^^] 

'  (5.19) 

-f  {[2  (a.6.6,  ^  u.B.B^  .  a^B.B.)  -  3  ^,^^\](F^^   F5) 

-  2  B.B.g,  F,,  --  (a.  .6,  +  a.,  6.  +  a.,6.)F,}  , 

1  J  k   4   a   ij  k    Ik  J    jk  1   6' 

'When  this  expression  is  substituted  into  the  cubic  form,  algebraic  simplification  yields 

2  2  ^ 

?      c    9^  In  f.  c  N  Z-;*     _  _         _ 

N  RT    I    ^-g-i  Z.7  Z^  =  RT  [  -  I   -^  .  3  Z  (6  F^)^  .  2(6  F^)^] 
i,j,k=1   J   k     -^  1  =  1   n. 

^  (5.20) 
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+  ^  [(6a  B^  -  3B^)(F2  +  F^)  -  2^'^?^    -  3  a  6  F^]  , 


where 


Z  s  ;_  Z.  ,  (5.21) 

i 


a  =   I   a.Z.  ,  (5.22) 

i 


=  ^  B.  Z.  ,  (5.23) 

i 


a  s  -  A  y  a. .  Z.  Z.  .  (5.24) 

a  ^  j   ij   1   J 

(There  is  a  typographic  error  in  the  last  term  on  the  right  side  of  eq   (18)  in  the  paper  by 
Michelsen  and  Heidemann  [7].)  To  compute  dT/dV  along  the  stability  limit,  we  need  the  partial 
derivatives  of  qy  with  respect  to  T  and  V;  these  are  given  by 

3  q  3F  8F 

N-^.  RT  [(6.  .   ,.)   —.2    B.Bj  F,  ^] 

(5.25) 

9F    a.  .  3F  8F, 

+  |-  [b.B.^ —  ^     "■    (B.6.  -a. 6.  -a.BJxTT-]  , 

b'-ij3K     adK       ij    ij    J13K-' 

where  the  derivatives  of  the  Fj^  with  respect  to  the  dimensionless  volume  K  are  given  at  the 
end  of  this  section,  and 

3   q.  .                 N    6.  .  „      B.8. 

N  -.7:;^^  =   R    i   ^  +    (6.    +    B.)    F,    +    B.B.F,     1    +  4|  -^   (F,   +   F,) 


3T  ^      n.  '"^i        ^^j'      1         '^i'^j    1    ^        3T        b      '    3  6' 

3   a.  ,   F^           3(aa.)    B  .        3(aa.)    B. 
JJ  _5  _   r i_  _A  +  J_  _il   F 


(5.26) 


9Tb  '■3Tb  3T  b    ■*      6 


where 


3   a..  (R   T    .)^  n.      1/2  -   m. 

^^  ^        P   ^  '  ^ci  (T  T    .)^^^ 

ci  ci 


9  a,.        1/2   C.  •  3   a..  9   a.. 

M '-^-TT^    [^^  a,,    .    a,,    ^r=^]  (5.28) 


3T  ^        1/2    ^    3T  jj        "ii      9T 

ii    jj 


9    (aa   )  c  9a.. 

-9f^  =   N    .^,    "j   -9T^  (5.29) 
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^ '  ?  .,i,  "'"j  ^  •  ''•'°' 


The  quadratic  forms  required  to  evaluate  dT/dV  can  be  simplified  to: 

NZ^^Z  =  RT  (2iZF^^  W  B^p^  F^^ ) 

^f  t^'  ^K3  -^^K5  ^  ^^'  ~2  a  B)  F^^^: 
and 

T  ^?  A  -  -  -2    2 

NZ-^Z  =  rUI  —  -^SBZF^  +  B   F^] 

i   i 


b  t^T  ^   (^3  '   h^    -  ^T  ^  -  2  "t  ^  ^6^  • 


;5.3i) 


(5.32) 


where 

8F. 

^T  "  If  '  ^^-^^^ 

c   5(aa.) 

a„  =  I   .T.  ^  Z   ,  (5.35) 

T    .^^   3T     1 

c   „ 
^T^i,l=1  -ir^^i^J   •  ^5.36) 

Dividing  (5.31)  by  (5.32)  gives  dT/dK  on  the  stability  limit  which  is  different  from  dT/dV  by 
a  factor  of  Nb. 

To  obtain  the  derivative  of  the  cubic  form  with  respect  to  K  along  the  stability  limit, 
differentiating  C  (5.20)  directly  gives 

2  3 

?  dr  c  N  Zf     _ 

N^  ^  =  RT^  [  -   I  —^   +  3  Z  (6  F^)^  +  2  (6  F^)^] 

i=1   n . 

1 

2   2 
c   "   Z   Z      _        _  __ 

■^  3RT  [  -  I     2 ^  ^  ^^V   ^  ^^^  ^K^^l  ^  ^^^   ""l^Kl 

i=l    n 

+  2  6^  6^  F^  ^  28^  F^  F^^  ] 

+  ^   1(6  aai^  -  3a8^)(F^2  *  ^k6^  ^  [^((aa)^^  2^  +  2  aaSBj^]         (5.37) 
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where 


-3   ia^^^   +   3a6^6j^)]    (F^   +   F^) 


-2    [(aj^e^   +   3ae^    Bj^)F^   +   ag^   F^^^; 


-3    [l(aa)^6   +  ^^^K^^e   *   ^^^   ^K6 ^ ^    ' 


T,  =-f.  (5.38) 


dZ.  ^^  dZ 

^Ki^     diT  '    ^        component   of  —   ,  (5.39) 


Z,   -f  =    j     Z^.    .  (5.n0) 

1  =  1 


1  =  1 


^  da  _  ^  dT  re;   ils>^ 

^K        dK        3T   dK    '  ^b.^^; 


1=1 


d(aa.)        3(aa.)     ,_, 
}^_  _  i_  oT 

dK        "        9T        dK    ' 


da..        0   a.  .     ._ 
y.  ^  iJ    dT 

dK  9T       dK 


(5.iJ'4; 


(5.46) 


When  tracing  critical  lines  in  regions  where  no  data  are  available  to  indicate  whether  or 
not  critical  points  actually  exist,  it  is  useful  to  be  able  to  compute  the  mathematical 
condition  for  stability,  i.e.,  that  the  quartic  form  is  positive.   The  quartic  form  for  the 
TCC  EOS  was  evaluated  by  taking  the  directional  derivative  of  the  cubic  form  along  the 
critical  displacement  vector;  that  is,  apply  the  operator 

^  "   ?   Z.  ,4-  (5.47) 

dn    .  ,   1  dn . 
1  =  1      1 
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to  eq      (5.20).      The  result    is 


3      4 
c     N^  Z 


i=1      r\i 

1 

+  I   [    12    (a6^   -   2   ai^  +    S^)(?2  +   F^)                                                 (5.48) 

-  8    (2   aB^   -    6^)   Fj^  +  5   6^   F^]    . 

5.1   Functions  of  the  Dimensionless  Volume  K  and  Their  First  Order  Derivatives 

Let  'A  =  6i  -  52,  then 

F^  =  1/(K-1)  , 

2   A  "-K+e^  K+5^''  ' 

R   =  1  r  f_J_i^  -f ^-1^1 

3   A  '-  '■K+S  ''  ''K+5  ''  ^    ' 

F  -  1  r  r-!L_i3  -f  ^2  >3i 

i)    A  '■  ''K+5  ■'  '■K+e  ''  J  ' 

9   A  '-  '■K+5  ''  ^K+S^    J  ' 
F^  =  4  In  '^ 


5   A    'K+&^' 


^6  =  ^2  -  ^5  • 


dF. 

and  their  first  order  derivatives  -rrr-  =   F,,.  are 

dK     Ki 


F   =  -  F  2 
KI      "l   ' 


K2   A     (K+6^)^    (K+S^) 


2  2 

F    =  1  [  -  1 . ?-^]    , 

K3   A     (K+6^)^  {K*6^)^ 
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X4 


H  4- 


K5     (K-H6^)(K+5^) 


F  ,  =  F 
K5    K2 


:<5  • 


5.   How  to  Use  the  Critical  Point  Software 

Two  subroutines  have  Deen  written  to  compute  critical  points.  CUBCR3  uses  a 
quasi-Newton  volume  iteration  and  CUBCRPT  uses  a  Newton-Haphson  volume  iteration  in  which 
the  volume  derivative  of  the  cubic  for:^  is  evaluated  analytically  (eq  5.37).  The  calling 
arguments  are  the  same  for  either  of  these  subroutines,  i.e., 

CALL  CUBCRPT  (NTYPE,  NC,  MXC,  X!i,  T,  ?,  V,  IPR,  CONV), 


On  entry, 

NTYPE 


NC 
MXC(I) 


XN(I) 
T 


IPR 

On  return, 
T 
P 
V 
CONV 


specifies  type  of  TCC  EOS.   Current  options  are: 

1  for  SHK  EDS  (Soave  [13]) 

2  for  PR  £03  (Peng  and  Robinson  [10]) 
#  of  mixture  components 

identification  number  for  each  mixture  component.   The  identification  number 

corresponds  the  position  of  the  component's  data  in  the  arrays  in  subroutine 

CUBDAT.  Check  a  source  listing  to  determine  the  pure  components  currently 

available. 

mole  numbers  (n,  n2,  ...,  n^). 

initial  guess  of  critical  temperature  (K).   If  set  to  0.0  an  initial  guess  is 

made  internally. 

initial  guess  of  critical  volume  (em3),  if  set  to  0.0  an  initial  guess  is  made 

internally. 

print  flag  for  debugging.  Set  to  0  for  no  output.   Set  to  1  for  summary  of 

iterations  written  to  unit  5. 

critical  temperature  (K) 
critical  pressure  (bar) 
critical  volume  (cm3) 
convergence  flag  (logical). 
=  .TRUE,  if  convergence  achieved 
=  .FALSE,  if  iteration  fails 


Externals, 


CUBPROP:   CUBDAT,  FVOL,  STABLIM,  CUBMXT ,  ZSUM,  CFORM,  S: 
CUBMX,  CUBP. 
LAS:   SUMVC,  DOT 
FORTRAN:   DABS,  DSIGN 
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:k,  zksum,  dck. 


Usage:        The  variables  XN,  T,  P,  and  V  are  double  precision.   The  subroutines  which  must 
be  linked  to  the  driver  are  located  in  the  following  files:   CUBCRIT,  CUBPROP, 
BIP,  DSPMAT,  and  LAS. 

The  software  is  structured  so  that  all  the  pure  component  data  that  are  required  to 
compute  a  mixture  critical  point  are  contained  in  subroutine  CUBDAT,   This  scheme  allows  the 
user  to  simply  specify  a  set  of  identification  numbers  for  the  mixture  components,  and  the 
number  of  moles,  of  each  component  in  order  to  compute  a  critical  point.  However,  to 
determine  the  critical  point  of  a  mixture  containing  a  component  (or  components)  for  which  the 
pure  component  data  are  not  contained  in  CUBDAT,  one  must  modify  the  software  to  include  the 
required  data.   The  pure  component  data  which  are  required  for  each  mixture  component  are:  a 
critical  temperature  (K) ,  critical  pressure  (bar),  and  the  acentric  factor.   The  arrays  in 
CUBDAT  should  be  dimensioned  as  follows: 

PARAMETER  (N=#  of  pure  components  whose  data  are  stored  in  CUBDAT,  NX=N*(N+1 )/2) 

DOUBLE  PRECISION  PC(N) ,  TTC(N) ,  W(N) ,  XI(NX) 

In  addition  to  the  pure  component  data,  any  binary  interaction  coefficients  (Cij)  which 
are  not  1.0  are  set  in  the  subroutines  SRKBIP  (for  the  SRK  EOS)  or  PRBIP  (for  the  PR  EOS). 
Just  follow  the  example  of  the  values  already  set  there. 

Storage  for  all  arrays  which  may  be  variably  dimensioned,  depending  on  the  number  of 
mixture  components,  is  set  in  subroutine  CUBCRPT  or  CUBCR3.   The  proper  dimensions  are  the 
following: 

PARA^4ETER  (N=#  of  mixture  components,  N1=N+1,  NX=N*  N1/2) 

INTEGER  KPVT(N) 

DOUBLE  PRECISION  AIJ(NX),  AIJT(NX),  AIT(N),  ALI(N),  ATIJ(NX),  AW(N),  BI(N),  BEI(N), 

Q(NX),  QA(N1,N),  TCA(N),  VCA(N),  Z(N),  ZK(N1),  ZZ(N) 

6.1   Verification  of  the  Software 

The  accuracy  of  the  software  was  verified  by  a  twofold  procedure: 

1)  the  results  of  a  binary  critical  point  calculation  were  compared  to  results  obtained 
from  an  independent  code  which  solved  the  traditional  critical  point  equations  of 
Gibbs; 

2)  the  software  was  checked  for  internal  consistency  by  comparing  selected  values  of 
derivatives  which  are  computed  using  the  formulas  given  in  this  note,  with  the  values 
obtained  using  different  formulas. 

Calculation  of  binary  critical  points  was  done  for  the  nCi5-C02  mixture  using  the  Peng 
and  Robinson  equation  of  state.   A  detailed  consideration  of  the  phase  diagram  for  nCi5-C02  as 
predicted  by  the  PR  EOS,  along  with  formulas  for  calculating  critical  points,  can  be  found  in 
Hong  [16]. 

Figure  6.1  shows  the  result  of  a  complete  tracing  of  the  critical  line  using  a  binary 
interaction  coefficient  <;  =  0.919,  and  the  following  pure  component  data: 

nCi6  CO2 

Tc(K)            717.0  30^.21 

Pc(MPa)             l.iJ2  7.38M 

w                  0.746  0.225 
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Figure  6.1   Calculated  critical  locus  for  the  hexadecane  +  carbon  dioxide  system 
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The  gas  constant  was  set  to  0.00831'J3'<  MPa'L/(mol -K) .   A  benchmark  result  for  x(nC-|5)  =  0.99 
follows: 

Tc  =  716.701292254  K 

Vc  =  1 .27759526809  L/raol 

Pc     =     1  .i46980319713  MPa 

The  independent  calculations  both  produced  the  above  result  to  the  accuracy  presented. 

The  following  consistency  checks  were  performed  on  multi-component  mixtures  containing  up 
to  12  components: 

1)  Values  of  the  cubic  form  computed  using  eq  (4.10)  were  compared  to  values  obtained 
from  the  "numerical"  version,  eq   (4.11); 

2)  Values  of  dT/dV  computed  using  eq  (4,18)  were  compared  to  values  obtained  from  a 
difference  formula; 

3)  Values  of  dZ/dV  obtained  by  solving  the  system  of  eq  (4.24)  were  compared  to  values 
obtained  from  a  difference  formula; 

4)  Values  of  dC/dV  computed  using  eq  (4.20)  were  compared  to  values  obtained  from  a 
difference  formula. 

Generally,  6  to  8  significant  figures  of  agreement  were  obtained  from  these  comparisons. 

7.    Summary  of  FORTRAN  Subprograms 

1 )   CUBDAT  computes  the  temperature  independent  part  of  the  pure  component  parameters  and 

mixture  cross  terms.   The  calling  sequence  is 

CALL  CUBDAT  (NTYPE,  NC,  MXC ,  G,  TCA,  VCA,  AIJ,  AW,  HI,  R). 

On  entry, 

NTYPE     specifies  type  of  TCC- EOS.   Current  options  are:   1  for  an  SRK  EOS,  2  for  a  PR 

EOS. 
NC        #  of  mixture  components. 
MXC(I)    identification  #  of  each  mixture  component. 

On  return, 

G(I)      G(1)  =  6i,  G(2)  =  62 

TCA(I)    critical  temperatures  (K)  of  pure  mixture  components 

VCA(I)    critical  molar  volumes  (cm3/mol)  of  pure  mixture  components 

AIJ(IP)   upper  triangular  part  of  the  symmetric  matrix  of  a^^j  packed  by  columns,  i.e., 

IP  =  i  +  j(j-1)/2.   In  unpacked  notation, 

AlJd.I)  =  fiaR^T^i/Pci 

AIJ(I,J)  =  Cij[AIJ(I,I)  »  AIJ(J,J)]l/2 
AW(I)     =  m(a)i) 
Bid)     =  %   RTci/Pci 
R        gas  constant  =  0.00831434  MPa'L/(mol -K) 

Externals, 

BIP:   3RKBIP,  PRBIP 
FORTRAN:   DSQRT 
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Usage:   CUBDAT  is  only  called  once  for  any  set  of  critical  point  or  stability  limit 
calculations  on  a  given  mixture  since  it  computes  quantities  which  are  independent  of  (T,  V, 
rii). 

2)   CUBMX  computes  the  temperature  and  composition  dependent  coefficients.   The  calling 
sequence  is 

CALL  CUBMX  (NC,  T,  XN,  XNT,  TCA,  AIJ,  AW,  BI,  ATIJ,  ALI,  BEI,  A,B). 

On  entry, 

NC  #  of  mixture  components  (c) 

T  temperature  (T)  in  kelvins 

XN(I)  mole  numbers  (n-|  ,  n2,  ...,  n,-,) 

c 
XNT       N  =   y  n. 
i  =  1  ' 
TCA(I)     Tci 

AIJ(IP)   temperature  independent  part  of  ay  (=a!.) 

AW(I)     m  (wi) 
BI       bi 

On  return, 

ATIJ(IP)   aj^j  in  packed  storage  of  upper  triangular  part 

ALKI)  a--    =  TT-     I     n.    a.  . 

1        Na    .^^      J      ij 


BEKI)  6i   =   bi/b 


1  ? 

a  =  — r       2,       n.n  .a.  . 

N^    i,j=1       '    J    ^J 


1  ? 

B  b   =  ^  I        n.b. 

Externals, 

FORTRAN:   DSQRT 
Usage:   CUBDAT  must  be  called  before  CUBMX. 

3)   CUBMXT  computes  the  temperature  derivatives  of  the  EOS  coefficients.   The  calling 
sequence  is 

CALL  CUBMXT  (NC,  T,  XN ,  XNT,  TCA,  AIJ,  AW,  ATIJ,  AT,  AIT,  AIJT). 

On  entry, 

NC        c 
T        T 
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XN(I)  rii 

XNT  N 

TCA(I)  Tci 

AIJ(IP)  a'ij 

AW(I)  m(a)i) 

ATIJ(IP)  aij 


On  return, 


c  da.. 


d   (aa. )  c  da.. 

^"^^)        -If-^  =  Jf  .1  "j  -dT 


da..        a..  ,      da..  .      da.. 

dT  2  a .  .      dT  a .  .        dT 

11  JJ 


where 


da..  -.1/2  -m. 
11           t       Ii         fi         (    T    ■)        1 1  r  1 
=   a! .    {1 +m. [1    - 


dT  ^^  ^  T    .  (T  T    .)^^^ 

ci  ci 

upper  triangular  part  of  symmetric  matrix  is  stored  in  packed  form. 

Externals, 

FORTRAN :   DSQRT 
Usage:         CUBDAT  and  CUBMX  must  be  called  before  CUBMXT. 

4)   FVOL  computes  the  functions  of  the  dimensionless  volume  (V/Nb)  and  their  1st  order 
derivatives.   The  calling  sequence  is 
CALL  FVOL  (V,  XNT,  B,  G,  FV) 

On  entry, 

V  system  volume  (V)  in  cm3 

XNT  N 

B  b 

G(I)  G(l)  =  6i,  G(2)  =  62 

On  return, 

FV(I)      FV(1)  =  Fi,  FV(2)=F2,  FV(3)=F3,  FV(i|)=Fij,  FV(5)=F5,  FV(6)=F5,  FV(7)=Fki, 
FV(8)=Fk2,  FV(9)=Fk3,  FV(10)=FKij,  FV(11)=Fk5,  FV(12)=Fk6.  FV(13)=F9 

Externals, 

FORTRAN:   DLOG 

32 


Usage:   CUF3DAT  must  be  called  before  FVOL.   The  argument  B  may  be  supplied  either  by  a  call  to 
CUBMX  or  simply  by  summing: 

b  =  -  I  n.b. 
N  .  ,   1  1 
1  =  1 

5)   CUBP  is  a  double  precision  function  which  returns  the  pressure.   The  calling  sequence 
is 

P  s   CUBP  (T,  VM,  G,  A,  B,  R) 

On  entry, 

T  T  (K) 

VM  V/N  (cm3/mol) 

G(I)  6i ,  62 

A  a 

B  b 

R  gas  constant 

On  return, 

CUBP      system  pressure  (P)  in  bar 
Usage:   CUBDAT  and  CUBMX  must  be  called  before  CUBP. 


6)   CUBDPV  is  a  double  precision  function  v;hich  returns  TTT?  j-r  m 

DPV  =  CUBDPV  (T,  VM,  G,  A,  B,  R) . 

The  argument  list  is  identical  to  that  for  CUBP. 


7)   CUBF  computes  fugacity.   The  calling  sequence  is 
CALL  CUBF  (I,  T,  V,  XN ,  XNT,  G,  A,  B,  ALI,  BEI,  F,  R) 


On  entry, 

I 

index  for  mixture  component  i 

T 

T 

V 

V 

XN(I) 

"i 

XNT 

N 

G(I) 

61  ,  62 

A 

a 

B 

b 

ALI(I) 

ai 

BEI(I) 

6i 

R 

gas  constant  R 

On  return. 

F 

RT  in  fi 
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Externals, 

FORTRAN:   DLOG 
Usage:   CUBDAT  and  CUBMX  must  be  called  before  CUBF 

8)  QIJ  computes  matrix  element  Qij.   The  calling  sequence  is 
CALL  QIJ  (I,  J,  T,  FV,  XN,  XNT,  ATIJ,  A,  B,  ALl,  BEl,  Q,  R) 

On  entry, 

I  index  for  mixture  component  i 

J  index  for  mixture  component  j 

T  T 

FV(I)  Fi,  FKi 

XN(I)  n^ 

XNT  N 

ATIJ(IP)  aj^j  (packed  storage) 

A  a 

B  b 

ALI(I)  ai 

BEI(I)  6i 

R  gas  constant 

On  return, 

8  In  f. 
Q        RT  — ^ — - 
J 

Usage:   CUBDAT,  CUBMX,  and  FVOL  must  be  called  before  QIJ.   This  routine  is  called  by  PAKQ  to 
assemble  Q. 

9)  PAKQ  assembles  the  upper  triangular  part  of  the  symmetric  Q  matrix  into  packed  form, 
the  calling  sequence  is 

CALL  PAKQ  (NC,  T,  FV,  XN ,  XNT,  Q,  ATIJ,  A,  B,  ALI,  BEI,  R) 


On  entry, 

NC 

#  of  mixture  components 

T 

T 

FV(I) 

Fi  and  Fki 

XIJ(I) 

"i 

XNT 

N 

ATIJ(IP)  aij 

A  a 

B  b 

ALIO  ai 

3EI(I)  6i 

R  gas  constant 

3>i 


On  return, 
Q(IP) 


packed  ^ 


Externals, 

CUBPROP:   QIJ 
Usage:   CUBDAT,  CUBMX,  and  FVOL  must  be  called  before  PAKQ.   PAKQ  is  called  by  DETQ  which 
computes  A,  and  by  DZK  which  computes  Zy^. 

10)   DETQ  assembles  g   by  calling  PAKQ,  then  computes  its  determinant  using  the  LINPACK 
routines  DSPFA  and  DSPDI.   A  flag  can  be  set  to  indicate  that  the  null  vector  is  to  be 
computed,  in  which  case  DSPCO  is  called  instead  of  DSPFA.   The  null  vector  is  normalized  to 
unit  length.   The  calling  sequence  is 
CALL  DETQ  (NC,  T,  FV ,  XN ,  XNT,  Q,  KPVT,  Z,  INL,  RCOND,  DTQ,  ATIJ,  A,  B,  ALI,  BEI,  R) 


On  entry. 

NC 

//  of  mixture  components 

T 

T 

FV(I) 

Fi  and  FKi 

XN(I) 

"i 

XNT 

N 

INL 

null  vector  computation 

ATIJ(IP) 

A 

B 

ALI(l) 

BEI(I) 

R 


INL  =  0  for  no  null  vector  computation 

INL  =  1  for  null  vector  computation 

When  INL  is  set  to  1 ,  it  is  assumed  that  the  £  matrix  has  already  been 

determined  to  be  singular. 


'ij 


gas  constant 


On  return, 
Q(IP) 
KPVT(I) 
Z(I) 


RCOND 

DTQ 


information  from  the  factorization  of  £  by  DSPFA  or  DSPCO 

array  used  by  LINPACK  routines  to  store  pivot  information 

null  vector  components  Z^^.   Z(I)  is  computed  by  LINPACK  routine  DSPCO,  and 

contains  an  approximate  null  vector  when  £  is  nearly  singular,  i.e.,  when  the 

reciprocal  condition  number  is  on  the  order  of  machine  roundoff. 

reciprocal  condition  number  computed  by  DSPCO 

A=det  I  a  I  computed  by  DSPDI 


Externals, 

CUBPROP:   PAKQ 

LINPACK:       DSPCO,    DSPFA,    DSPDI 
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FORTRAN:   I  DINT,  DSQRT 
Usage:   CUBDAT,  CUBMX,  and  FVOL  must  be  called  before  DETQ.   This  routine  is  called  by  QNT  to 
find  T  on  the  stability  limit,  then  compute  !_. 

11)   QNT  performs  the  quasi-Newton  Iteration  on  temperature  to  solve  A(T)  s  det  |£|  =  0 
for  a  given  volume.   The  calling  sequence  is 

CALL  QNT  (NC,  T,  FV,  XN ,  XNT,  RDT,  TOLT,  ITMX,  IPR,  IT,  IFT,  Q,  KPVT,  Z,  TCA,  AIJ,  AW,  ATIJ, 
EI,  ALI,  BEI,  R) 


On  entry, 
NC 
T 

FV(I) 
XN(I) 
XNT 
RDT 
TOLT 

ITMX 
IPR 


TCA(I) 

AIJ(IP) 

AW(I) 

Bid) 

R 


#  of  mixture  components 
initial  guess  of  TCTq) 


N 

6T/T  where  6T  is  the  increment  used  to  compute  the  derivative  3A(TQ)/tJT 

iteration  convergence  tolerance 

I  Ti  +  i  -  Ti  |/Ti  <  TOLT 

maximum  number  of  iterations 

print  flag: 

0  for  no  printing 

1  for  printing  of  iteration  progress  on  unit  6 
Tci 

m(wi ) 

bi 

gas  constant 


On  return, 
T 
IT 

IFT 


temperature  on  stability  limit  (A=0) 
number  of  iterations  to  convergence 
iteration  flag: 

0  successful  iteration 

1  iteration  failed 
Q(IP)     factored  £  (from  DSPFA) 
KPVT(I)    pivots  (LINPACK  routines) 
Z(I)      not  used 

ATIJ(IP)   aj^j  at  T  on  stability  limit 
ALI(I)    tti  at  T  on  stability  limit 
BEKI)    Si 


Externals, 

CUBPROP:   CUBMX,  DETQ 
FORTRAN:   DABS,  DSIGN 
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Usage:   CUBDAT  and  FVOL  must  be  called  before  QNT.   CUBMX  and  DETQ  are  called  from  within  the 
iteration  loop  on  temperature  in  QNT.   This  routine  is  called  by  STABLIM. 

12)   STABLIM  locates  the  temperature  on  the  stability  limit  when  the  volume  is  fixed. 

Then  the  null  vector  is  computed.   On  the  1st  call  to  STABLIM  the  direction  of  Z^    is 

determined  by  requiring  Z-j  >  0.   On  subsequent  calls  to  STABLIM,  for  a  given  critical  point  or 

T 
stability  limit  calculation,  the  direction  of  Z_.  is  determined  by  requiring  Z._  Z.  >  0. 

The  calling  sequence  is 

CALL  STABLIM  (NC,  T,  FV,  XN ,  XNT,  RDT,  TOLT,  ITMX,  IPR,  IT,  IFT,  IFC,  Q,  KPVT,  Z,  ZZ,  TCA, 

AIJ,  AW,  ATIJ,  BI,  ALI,  BEI,  A,  B,  R) 


On  entry, 
NC 
T 

FV(I) 
XN(I) 
XNT 
RDT 
TOLT 

ITMX 

IPR 


IFC 


TCA (I) 

AIJ(IP) 

AW(I) 

BI(I) 

R 


#  of  mixture  components 

initial  guess  of  T  on  stability  limit 


N 

6T/T,  where  6T  is  the  increment  used  to  compute  the  derivative  3A(Tq)/8T 

iteration  convergence  tolerance,  i.e., 

I  Ti  +  i  -  Til  /Ti  <  TOLT 

maximum  number  of  iterations 

iteration  print  flag: 

0  for  no  printing 

1  for  printing  of  iteration  progress  on  unit  6 
first  call  flag 

1  first  call;  set  Zj  >  0 

T 
0  not  first  call;  require  Z.  ,  Z.  >  0 
^     -1-1  -1 

Tci 

m((i)j^ ) 

bi 

gas  constant 


On  return, 
T 

IT 
IFT 


Q(IP) 

KPVT(I) 

Z(I) 

ZZ(I) 

ATIJ(IP) 

ALI(I) 


temperature  on  stability  limit 

#  of  iterations  for  convergence 

iteration  convergence  flag--not  currently  used--if  iteration  fails  in  QNT, 

execution  is  stopped  in  STABLIM. 

factored  g   (from  DSPCO) 

pivot  information  (from  DSPCO) 

li 

Z^j^-I  (from  previous  iteration  on  volume) 
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BEKI)  Bi 

A  a 

B  b 

Externals, 

CUBPROP:  QNT,  CUBMX,  DETQ 

FORTRAN:  DABS 

Usage:   CUBDAT  and  FVOL  must  be  called  before  STABLIM. 

13)   ZSUM  computes  the  summations  involving  the  null  vector  components  which  are  used  to 
compute  the  quadratic,  cubic,  and  quartic  forms.   The  calling  sequence  is 

CALL  ZSUM  (NC,  XN,  XNT,  Z,  ATIJ,  ALI,  BEI,  A,  AIT,  AIJT,  ZS) 

On  entry, 

NC  #  of  mixture  components  (c) 

XN(I)  ni 

XNT  N 

Z(I)  null  vector  components  (Z-|  ,  Z2 1q) 

ATIJ(IP)  aj^j  (packed  storage) 

ALI(I)  tti 

BEI(I)  6i 

A  a 

AIT(I)  d(aai)/dT 

AIJT(IP)  d(aij)/dT  (packed  storage) 


On  return. 


2 
c  Z.  n. 


ZS(I)     ZS(1)  =  I     ^^  wher-    "   ^ 


..,.^,  e  y    ^ 
1=1  ■'1 


c  Z.3 


ZS(2)  =  I     -ir 


c   Z.'* 
ZS(9)  =  I     ^ 


c 
zs(3)  =  z  =  I    z. 
i=1   ^ 


c 
ZS(4)  =6=1   B.Z 
i=1 
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ZS(5)  =  a  =  I     a  Z 
i  =  1 


-   1    % 
ZS(6)  =  a  =  -   y   a, ,  Z.Z. 


c   d  (aa. ) 
ZS(7)  =  a,  =  I     — Ht^Z. 
1=1 


ZS(8)  =l^=        l^   Z.Z 
i.J  =  l 

Usage:   CUBDAT,  FVOL,  STABLIM,  and  CUBMXT  must  be  called  before  ZSUM. 

1^)   QFORM  computes  the  quadratic  form.   The  calling  sequence  is 
CALL  QFORM  (XNT,  R,  T,  A,  B,  ZS,  FV,  QF) 

On  entry, 

XNT  N 

R  gas  constant 

T  temperature  on  stability  limit 

A  a 

B  b 

ZS(I)  sums  involving  Z 

FV(I)  Fi  and  F^^ 

On  return, 

QF        zT  g  Z 
Usage:   CUBDAT,  FVOL,  STABLIM,  CUBMXT,  and  ZSUM  must  be  called  before  QFORM. 

15)  CFORM  computes  the  cubic  form.   The  calling  sequence  is 
CALL  CFORM  (XNT,  R,  T,  A,  B,  ZS,  FV,  CF) 

On  entry, 

Same  as  QFORM 

On  return, 

?     a^A 

OF  I  .      I   ^,     Z.Z.Z, 

.  .  ,  ,  3n. 3n . 3n,    i  J  k 

Usage:   same  as  QFORM 

16)  DFORM  computes  the  quartic  form.   The  calling  sequence  is 
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CALL  DFORM  (XNT,  R,  T,  A,  B,  ZS,  FV,  DF) 
On  entry,  same  as  QFORM. 

On  return, 

c  4 

DF  I  i: — ,  ^/  .  -  Z.Z.Z,  Z^ 

.  .  ,  -,  ,  9n.  3n  .dn,  dn,   i  j  k  1 
i,j,k,l=l    1  J   k  1 

Usage:   same  as  QFORM. 

17)   SLDTK  computes  dT/dK  along  the  stability  limit.   The  calling  sequence  is 
CALL  SLDTK  (R,  T,  A,  B,  AT,  ZS,  FV,  DTK) 


On  entry, 

R 

gas  constant 

T 

T 

A 

a 

B 

b 

AT 

da/dT 

ZS(I) 

sums  involving  Z 

FV(I) 

Fi  and  Fki 

On  return, 

DTK       dT/dK 
Usage:   CUBDAT,  FVOL,  STABLIM,  CUBMXT,  and  ZSUM  are  called  before  SLDTK. 

18)   QKMAT  assembles  the  upper  triangular  part  of  the  symmetric  matrix  d£/dK  in  packed 

form.   The  calling  sequence  is 

CALL  QKMAT  (NC,  T,  FV,  XN ,  XNT,  QK ,  ATIJ,  A,  B,  ALI,  BEI,  AT,  AIT,  AIJT,  DTK,  R) 

On  entry, 

NC  #  of  mixture   components 

T  temperature  on  stability   limit 

FV(1) 

XN(I)  n^ 

Xt>lT  N 

ATIJ(IP)      a^j 

A  a 


B 

b 

ALI(I) 

ai 

BEKI) 

h 

AT 

da/dT 

AIT(I) 

d(aai)/dT 

AIJT(IP) 

d  Sij/dT 

DTK 

dT/dK 

i»0 


R        gas  constant 

On  return, 

QK(IP)    d£/dK  =  8g/3K  +  0£/8T) (dT/dK) 
Usage:   CUBDAT,  FVOL,  STABLIM,  CUBMXT,  ZSUM,  and  SLDTK  must  be  called  before  QKMAT 

19)   DZK  computes  dZ^/dK  along  the  stability  limit.   In  assembling  the  rectangular  matrix, 
the  bottom  row  (Z^"^)  is  scaled  so  that  the  largest  component  of  Z  is  the  same  magnitude  as  the 
largest  component  of  Q.      The  rectangular  system  is  solved  using  a  simple  Gaussian  elimination 
algorithm  with  partial  pivoting,  i.e,  row  exchanges  only.   Thus,  at  some  point  in  the 
elimination  sequence,  a  zero  pivot  is  encountered  (because  £  is  singular)  and  the  bottom  row 
will  be  exchanged  with  the  row  containing  that  pivot.   The  elimination  sequence  is  completed, 
and  the  back  substitution  proceeds  as  if  the  system  were  square.   The  calling  sequence  is 
CALL  DZK  (NC,  T,  FV,  XN,  XNT,  Q,  ATIJ,  A,  B,  ALI,  BEI,  AT,  AIT,  AIJT,  DTK,  NCI,  MQA,  QA,  Z, 
ZK,  R) 


On  entry, 
NC 
T 

FV(I) 
XK(I) 
XNT 

ATIJ(IP) 
A 
B 

ALI(I) 
BEI(I) 
AT 

AIT(I) 
AIJT(IP) 
DTK 
NCI 
MQA 


IJ 


#   of  mixture  components  (c) 

temperature  on  stability  limit 

F 

n 

N 

a 

a 

b 

da/dT 
d(aui)/dT 
d  a^j/dT 


dT/dK 

c  +  1 

1st  dimension  of  QA  array 
QA(I,J)   work  space  used  for  rectangular  matrix.  Must  be  dimensioned  at  least  QA  (NCI, 

NC) 
Z(I)      Zi 
R        gas  constant 


On  return, 
Q(IP) 
ZK(I) 


Externals, 

CUBPROP 
LAS: 


PAKQ,  QKMAT 
MSAX,  RMSOL 
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FORTRAN:   DABS 
Usage:   CUBDAT,  FVOL,  STABLIM,  CUBMXT,  ZSUM,  and  SLDTK  must  be  called  before  DZK. 

20)   ZKSUM  computes  the  sums  involving  Zj^j^ .   The  calling  sequence  is 
CALL  ZKSUM  (NC,  XN ,  XNT,  Z,  ZK,  ATIJ ,  A,  ALI,  BEI,  ZS,  DTK,  ZKS) 

On  entry, 

NC  #  of  mixture  components 

XN(I)  ni 

XNT  N 

Z(I)  Zi 

ZK(I)  Zki 

ATIJ(IP)  aij 

A  a 

ALI(I)  tti 

BEI(I)  6i 

ZS(I)  sums  involving  Zj^ 

DTK  dT/dK 

On  return, 

c   Z^  Z  n. 

ZKS(I)    ZKS(l)  =  I     -^ — -   where  y  =  — 

i=1   y.  i   N 

1 


ZKS(2)  =  Z^   .  I      Z^. 
1=1 


ZKS(3)  =  B,  -     l^    B.Z^. 


c  8(aa  ) 

ZKS(^)  =  (aa)^=-   I   [aa.  Z^^.  .  Z.  -3^  - 
1=1 


ZKS(5)  =  (aa)j^ 

Usage:   CUBDAT,  FVOL,  STABLIM,  CUBMXT,  ZSUM,  SLDTK,  and  DZK  must  be  called  before  ZKSUM. 

21)  DCK  computes  the  derivative  of  the  cubic  form  along  the  stability  limit.   The  callinf 

sequence  is 

CALL  DCK  (XNT,  R,  T,  DTK,  A,  B,  AT,  ZS,  FV,  ZKS,  CK) 
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On  entry, 

XNT 

N 

R 

gas  constant 

T 

temperature  on  stability  limit 

DTK 

dT/dK 

A 

a 

B 

b 

AT 

da/dl 

ZS(I) 

sums  involving  Z 

ZKS(l)    sums  involving  Z^ 

On  return, 

CK        dC/dK 
Usage:   CUBDAT,  FVOL,  STABLIM,  CUBMXT,  ZSUM,  SLDTK,  DZK,  and  ZKSUM  must  be  called  before  DCK. 

22)  MSAX  computes  the  product  of  a  symmetric  matrix  and  a  vector.   The  calling  sequence 
is 

CALL  MSAX  (N,  A,  X,  B) 

On  entry, 

N        order  of  the  matrix 

A(IP)     upper  triangular  part  of  A  stored  in  packed  form,  i.e.,  A( J*( J-1 )/2+I)  =  A(I,J) 

X(I)      vector  to  be  multiplied  by  A 

On  return, 

8(1)      result  A  X  =  B 

23)  SUMVC  sums  the  components  of  a  vector.   The  calling  sequence  is 
CALL  SUMVC  (N,X,S) 

On  entry, 

N        size  of  vector 

X(I)      vector  to  be  summed,  x 

On  return, 

N 
S         y  X. 
i  =  l  ' 

24)  DOT  computes  the  dot  product  of  two  vectors.  The  calling  sequence 
CALL  DOT  (N,X,Y,S) 

On  entry, 

N  size   of   vectors 

>i3 


On 

entry, 

NA 

M 

N 

Ad.j) 

B(I) 

On 

return, 

B(I) 

X(I),Y(I)  vectors  to  be  dotted;  X.  X 

On  return, 

N 
S  y  X.Y. 

i  =  l  '    ' 

25)   RMSOL  solves  the  rectangular  system  A  2^  ^  M.'  "here  A  is  an  M  x  n  matrix.   It  is 
assumed  that  M  =  N+1  and  that  a  unique  solution  exists,  i.e.,  the  rank  of  A  is  N.   The  matrix 
is  factored  by  Gaussian  elimination  with  partial  pivoting.   The  calling  sequence  is 
CALL  RMSOL  (NA,  M,  N,  A,  B) 


1st  dimension  of  the  array  A(I,J) 

#  of  rows  in  ^  (assumed  to  equal  N+1) 

#  of  columns  in  A 
matrix  to  be  factored 
right  hand  side  B 


solution  X 

8.    References 

[I]  Courant,  R. ;  John,  F.  Introduction  to  calculus  and  analysis,  Vol.  2  New  York:  John  Wiley 
&  Sons;  197^1.  187. 

[2]   Dongarra,  J.;  Bunch,  J.  R. ;  Moler,  C.  B. ;  Stewart,  G.  W.  LIMPACK  user's  guide. 

Philadelphia,  PA:  SIAM  Publications;  1978, 
[3]   Gibbs,  J.  W.  On  the  equilibrium  of  heterogeneous  substances  (Oct.  1876-May  1877), 

collected  works,  Vol.  1.  New  Haven  CT:  Yale  Univ.  Press;  1928.  55. 
[4]   Heidemann,  R.  A.;  Khalil,  A.  M.  The  calculation  of  critical  points.  AIChE  J.  26: 

769-779;  1980. 
[5]   Hildebrand,  F.  B.  Methods  of  applied  mathematics,  2nd  ed.  Englewood  Cliffs,  NJ : 

Prentice-Hall  Inc.;  1965.  52. 
[6]   Michelsen,  M.  L.  Calculation  of  phase  envelopes  and  critical  points  for  multicomponent 

mixtures.  Fluid  Phase  Equilibria  M:  1-10;  1980. 
[7]   Michelsen,  M.  L.;  Heidemann,  R.  A.  Calculations  of  critical  points  from  cubic 

two-constant  equations  of  state.  AIChE  J.  27:  521-523;  1981. 
[8]   Modell,  M. ;  Reid,  R.  C.  Thermodynamics  and  its  applications.  Englewood  Cliffs,  NJ : 

Prentice-Hall  Inc.;  197^. 
[9]   Peng,  D.-Y.;  Robinson,  D.  B.  A  new  two-constant  equation  of  state.  Ind.  Eng.  Chem. 

Fundam.  15:  59-6iJ;  1976. 
[10]  Peng,  D.-Y.;  Robinson,  D.  B.  A  rigorous  method  for  predicting  the  critical  properties  of 

multicomponent  systems  from  an  equation  of  state.  AIChE  J.  23:  137-1^^1;  1977. 

[II]  Peters,  G. ;  Wilkinson,  J.  H.  Inverse  iteration,  ill-conditioned  equations  and  Newton's 
method,  SIAM  Review  21:  339-360;  1979, 


[12]  Prausnitz,  J.  M.  Molecular  thermodynamics  of  fluid-phase  equilibria.  Englewood  Cliffs, 

NJ:  Prentice-Hall  Inc.;  1969. 
[13]  Soave,  G.  Equilibrium  constants  from  a  modified  Redlich-Kwong  equation  of  state.  Chem 

Eng.  Sci.  27:  1197-1203;  1972. 
[14]  Tisza,  L.  Generalized  thermodynamics.  Cambridge,  MA:  The  M.I.T.  Press;  1977. 
[15]  Weinhold,  F.  Thermodynamics  and  geometry.  Physics  Today.  23;  March  1976. 
[16]  Hong,  G.T.;  Binary  phase  diagrams  from  a  cubic  equation  of  state.   Ph.D  Thesis.   MIT; 

1981  . 


^5 


NBS-lUA    (REV.    2-8C) 


U.S.    DEPT.    OF    COMM. 

BIBLIOGRAPHIC  DATA 

SHEET  (See  instructions) 


1.  PUBLICATION  OR 
REPORT  NO. 

NBS/TN-1313 


2.  Performing  Organ.  Report  No, 


3.  Publ  ication  Date 

March   1988 


4.  TITLE  AND  SUBTITLE 

ON  THE  CALCULATION  OF  CRITICAL  POINTS  BY  THE  METHOD  OF  HEIDEMANN  AND  KHALIL 


5.  AUTHOR(S) 

Brian   E.    Eaton 


6.  PERFORMING  ORGANIZATION  (tf  joint  or  other  thon  NBS,  see  in  struct  ion  sj 

NATIONAL  BUREAU  OF  STANDARDS 
DEPARTMENT  OF  COMMERCE 
WASHINGTON,  D.C.    20234 


7.  Contract/Grant  No. 


8.  Type  of  Report  &  Period  Covered 


9.  SPONSORING  ORGANIZATION  NAME  AND  COMPLETE  ADDRESS  (Street,  City.  State,  ZIP) 


10,  SUPPLEMENTARY  NOTES 


^      Document  describes  a  computer  program;  SF-185,  FlPS  Software  Summary,  is  attached. 


11.  ABSTRACT  (A  2Q0-word  or  less  factual  summary  of  most  significant  information.    If  document  includes  a  significant 
bi bliography  or  literature  survey,  mention  it  here) 

The  formulation  of  critical  point  criteria  by  Heidemann  and  Khalil  is  analyzed 
and  contrasted  with  the  original  formulation  of  Gibbs.   An  extension  to  the 
solution  technique  originally  used  by  Heidemann  and  Khalil,  and  later  improved 
by  Michelsen  and  Heidemann,  is  presented  along  with  its  detailed  implementation 
for  a  general  two-constant  cubic  equation  of  state.   Finally,  FORTRAN  software 
developed  for  these  computations  is  carefully  documented. 


12.  KEY  WORDS  (Six  to  twelve  entries;  alphabetical  order;  capitalize  only  proper  names;  and  separate  key  words  by  semicolons) 

algorithm;  critical  locus;  Helmholtz  free  energy;  mixtures;  Peng-Robinson  equation 
of  state;  thermodynamics 


13.  AVAILABILITY 

|_Xj  Unlimited 

I      I  For  Official  Distribution.    Do  Not  Release  to  NTIS 

fX]  Order  From  Superintendent  of  Documents,  U.S.  Government  Printing  Office,  Washington,  D.C. 
20402. 

'~      Order  From  National  Technical  Information  Service  (NTIS),  Springfield,  VA.    22161 


14.  NO.  OF 

PRINTED  PAGES 

52 


15.  Price 


USCOMM-DC    6043-P80 


w  US  GOVERNMENT  PRINTING  OFFICE:  1988— 575001.85.119  REGION  NO.  8 


NBS 


_  Technical  Publications 

Periodical 

Journal  of  Research — The  Journal  of  Research  of  the  National  Bureau  of  Standards  reports  NBS  research 
and  development  in  those  disciplines  of  the  physical  and  engineering  sciences  in  which  the  Bureau  is  aaive. 
These  include  physics,  chemistry,  engineering,  mathematics,  and  computer  sciences.  Papers  cover  a  broad 
range  of  subjects,  with  major  emphasis  on  measurement  methodology  and  the  basic  technology  underlying 
standardization.  Also  included  from  time  to  time  are  survey  articles  on  topics  closely  related  to  the  Bureau's 
technical  and  scientific  programs.  Issued  six  times  a  year. 

Nonperiodicals 

Monographs — Major  contributions  to  the  technical  literature  on  various  subjects  related  to  the  Bureau's  scien- 
tific and  technical  activities. 

Handbooks — Recommended  codes  of  engineering  and  industrial  practice  (including  safety  codes)  developed  in 
cooperation  with  interested  industries,  professional  organizations,  and  regulatory  bodies. 

Special  Publications — Include  proceedings  of  conferences  sponsored  by  NBS,  NBS  annual  reports,  and  other 
special  publications  appropriate  to  this  grouping  such  as  wall  charts,  pocket  cards,  and  bibliographies. 
Applied  Mathematics  Series — Mathematical  tables,  manuals,  and  studies  of  special  interest  to  physicisLs, 
engineers,  chemists,  biologists,  mathematicians,  computer  programmers,  and  others  engaged  in  scientific  and 
technical  work. 

National  Standard  Reference  Data  Series — Provides  quantitative  data  on  the  physical  and  chemical  properties 
of  materials,  compiled  from  the  world's  literature  and  critically  evaluated.  Developed  under  a  worldwide  pro- 
gram coordinated  by  NBS  under  the  authority  of  the  National  Standard  Data  Act  (Public  Law  90-3%). 
NOTE:  The  Journal  of  Physical  and  Chemical  Reference  Data  (JPCRD)  is  published  quarterly  for  NBS  by 
the  American  Chemical  Society  (ACS)  and  the  American  Institute  of  Physics  (AIP).  Subscriptions,  reprints, 
and  supplements  are  available  from  ACS,  1155  Sixteenth  St.,  NW,  Washington,  DC  20056. 

Building  Science  Series — Disseminates  technical  information  developed  at  the  Bureau  on  building  materials, 
components,  systems,  and  whole  structures.  The  series  presents  research  results,  test  methods,  and  perfor- 
mance criteria  related  to  the  structural  and  environmental  functions  and  the  durability  and  safety 
characteristics  of  building  elements  and  systems. 

Technical  Notes — Studies  or  reports  which  are  complete  in  themselves  but  restrictive  in  their  treatment  of  a 
subject.  Analogous  to  monographs  but  not  so  comprehensive  in  scope  or  definitive  in  treatment  of  the  subject 
area.  Often  serve  as  a  vehicle  for  final  reports  of  work  performed  at  NBS  under  the  sponsorship  of  other 
government  agencies. 

Voluntary  Product  Standards — Developed  under  procedures  published  by  the  Department  of  Commerce  in 
Part  10,  Title  15,  of  the  Code  of  Federal  Regulations.  The  standards  establish  nationally  recognized  re- 
quirements for  products,  and  provide  all  concerned  interests  with  a  basis  for  common  understanding  of  the 
characteristics  of  the  products.  NBS  administers  this  program  as  a  supplement  to  the  activities  of  the  private 
sector  standardizing  organizations: 

Consumer  Information  Series — Practical  information,  based  on  NBS  research  and  experience,  covering  areas 
of  interest  to  the  consumer.  Easily  understandable  language  and  illustrations  provide  useful  background 
knowledge  for  shopping  in  today's  technological  marketplace. 

Order  the  above  NBS  publications  from:  Superintendent  of  Documents,  Government  Printing  Office, 
Washington.  DC  20402. 

Order  the  following  NBS  publications — F/PS  and  NBSIR  's^rom  the  National  Technical  Information  Ser- 
vice, Springfield,  VA  22161. 

Federal  Information  Processing  Standards  Publications  (FIPS  PUB) — Publications  in  this  series  collectively 
constitute  the  Federal  Information  Processing  Standards  Register.  The  Register  serves  as  the  official  source  of 
information  in  the  Federal  Government  regarding  standards  issued  by  NBS  pursuant  to  the  Federal  Property 
and  Administrative  Services  Act  of  1949  as  amended.  Public  Law  89-306  (79  Stat.  1 127),  and  as  implemented 
by  Executive  Order  1 1717  (38  FR  12315,  dated  May  11,  1973)  and  Part  6  of  Title  15  CFR  (Code  of  Federal 
Regulations). 

NBS  Interagency  Reports  (NBSIR) — A  special  series  of  interim  or  final  reports  on  work  performed  by  NBS 
for  outside  spxjnsors  (both  government  and  non-government).  In  general,  initial  distribution  is  handled  hy  the 
sponsor;  public  distribution  is  by  the  National  Technical  Information  Service,  Springfield,  VA  22161,  in  paper 
copy  or  microfiche  form. 
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